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FLUID FLOW BETWEEN POROUS ROLLERS 


By G. I. TAYLOR and J. C. P. MILLER (Cambridge) 
[Received 5 May 1955] 
SUMMARY 

When a viscous fluid is entrained between two rollers which are separated by a 
small distance, a pressure is developed on the upstream side of the point of nearest 
approach and a suction on the downstream side. Theoretically these both become 
infinite when the rollers are in contact. If, however, the rollers are slightly porous 
these infinities are avoided. The analysis of this situation involves the solution of 
a differential equation which has some interest. The distribution of pressure is 
rather similar to that which occurs when impervious rollers are not in contact. 
The result may have industrial applications to paper-making machines and to the 
rollers used for painting walls. 

1. Introduction 

Iy certain industrial processes a viscous fluid is caught up between two 
rotating rollers or between a moving flat sheet and a roller which supports 
it. The pressure upstream of the point where the rollers most nearly touch 
first rises as this point is approached then falls to zero when the point is 
reached. Then after passing this point the pressure falls below that at 
great distances and subsequently rises. 

The distribution of pressure between cylindrical rollers which are not 
porous has been calculated (1) on the assumption that they are separated 
by a very narrow gap so that the equations o 
of the theory of lubrication can be used. The 
pressure p can then be regarded as a function 
of x only, where the origin of coordinates is 
chosen at the point midway between the centres 
of two equal cylinders of radius R, and the axis 
of z is the line joining their centres (Fig. 1). 





If U is the velocity of the surface of the 
rollers and 2h the small distance between them, 


the pressure p is 
4ul Rx 
be ae . ] 
P (27+ 2 Rh)? ") 


where p is the coefficient of viscosity. Since zis yg. 





Impervious rollers 
measured in the direction of motion the pressure rotating at the same speed 


: : ‘ ; R : in opposite directions. 
18 antisymmetrical, being a suction when vis — aes 


Positive, ie. downstream of the point where the gap is narrowest. This 
suction is a maximum at the point x = (3.Rh)! and is of magnitude 


(=P) nae = fy V6uU Rth-', 


(Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 2 (1956)] 
5002.34 K 
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If h is decreased (— p),,,,, increases till, when the rollers come into contact. 
the formula predicts an infinite valne. In fact (2) the suction increases til] 
at a certain value, depending on its physical properties, the fluid cavitates, 
There is however another possibility, 
The rollers themselves might be slightly 
porous and when the geometrical con- 
ditions are such that high suction would 
be produced by impermeable rollers, 
fluid might be sucked through their 
surfaces and so prevent the suction 
rising to such a level that cavitation 
would occur. It is of interest to discuss 
the flow in this case, partly to find out 
how porous the roller must be if cavita- 





tion is to be avoided but perhaps more 
Fic. 2. Porous rollers in contact. so because the analysis presents some 
interesting mathematical features. 
It will be assumed that the velocity W of the fluid through the porous 
surface of the rollers is related to the pressure p by the equation 
W kp -; (2) 
k is then a porosity coefficient. It has the dimensions of a length. The 


geometry of the situation is indicated in Fig. 2. 


2. Flow equations 

Since only the region close to the point of contact of the rollers will be 
considered it is sufficient to take 2t = 2?/R as the gap between them at 
the point x. 

According to the theory of lubrication the velocity of flow at the point 


(7,2) is 
- rw l “2 dp 
2u dx 
If U is the velocity of the surface of the rollers 
4 
A : dp (3) 
Su k2 dx 
and the total flow Q is 
t 
277 6 
Q 2 u dz i | x dp (4) 


R 12 Ru dx 
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act, The equation of continuity is 
til dQ 

- 2W 
ites, dz 
lity bret ; 
ht] so that, using (2), Q and W can be eliminated to give 
tly 


, ae are 

con | a (x6 =A 2kp 2al 0, (5) 
yuld L2u R3 da dx ph R 

lers und writing 

kR | 

net vy v(24k.h*)-*, y p ;)(24kR®) t (6) 
ul 

d ; dy ‘ 

sia 5) becomes dz, vy ae y+, (). (7) 


This equation may be transformed by the substitution 
| 

Dy2’ 

—s 

to a modified form of the Bessel function equation, viz. 

TOUS dy dy 

ede" ° a ' 


- The compl mentary function of (9) is 





The | 1’1;(€)+ B’L_s(€). (10) 


i 
lo complete the solution it is necessary to find a particular integral of (9). 


For this purpose consider Struve’s function (3, p. 329) 


————— 


L¢é\v+2m+1 
(2¢) 


ill be | mks pa P (m+ 3)P(m+-v 3)’ (11) 
. t which satisfies (3. ibid.) 

} \ : 4(3£)"+} (12) 
p r(s)Pw+3) 
I where Y written tor the operator 
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so that a particular integral of (9) is 
512-1 (4A) I(3) Ly (€)+ 2-23-10) (17) 

and the complete solution of (7) is 

y = (2€)[2-8P (3) PR) L,(€)+- 2-23-41 + AL (é)+ BL ,(6)]. (18) 
When the functions L;(&), J;(€), 1_y(£) are expanded in powers of é the only 
term which does not vanish when é + 0 is that which contains B. Since 
the physical conditions require that p > 0 as x + 00, or y = 0 when é = 0, 
therefore B = 0. We next consider whether any value can be chosen for 
A which will make y > 0 as >. Both J,() and L,(&) tend to infinity 
as € -> 00 but the difference between them is of a lower order of magnitude 
than either separately. This can be seen by writing down the expressions 
for L,(€) and /,(€) as definite integrals. These are (cf. 3, p. 329) 


hor 


24g =f 





L,(é) Po+hre) sinh(é cos @)sin?"6 dé (19) 
0 
and (3, p. 79) 
hor 
y(1l¢\v le 
Dé) neenn 7 | cosh(€é cos 6)sin?”6 dé; (20) 
0 
hence L,(é)—I,(€) = orate | e-€ 008 8 gin2vQ dQ (21) 


and the reason for the special character of this particular combination of 
L,(é) and [,(€) is obvious. 

When € is large only the part of the range of @ which lies close to }7 
contributes appreciably to the integral in (21). Writing ¢ = 47—6 this 
integral is = 

9 = | e-€sind cos?"d dd. (22) 

0 
To obtain the first two terms in the expansion of (22) in descending powers 
of € write sind = d—)¢d* = so that to the same order ¢ = 4+ }/3. Then 
(1—2v) 


1 
: ] 
J ~ | ¢ els +-(4—v)p?} dip ~ z+ 
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(23) 


Hence, from (21), 
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L,(é) and J,(€) occur in the combination of (21) and the complete solution is 
ni y = (2€)'f2-1F(4)P(3){ Ly (6) —4(6)} + 243414]. (25) 
The limiting form of (25) as € > co is found by substituting from (24) in 
i" 25). It is ) 9-1(¢4__ 3¢-i »—39-1¢) ; 
y ~ (26) 2"*3 iSo*—-aE 4)+2 13 é']. (26) 
only 
Since The term in €! vanishes leaving 
¥ y ~ (2€)8(2E)-= = (2€)-? = ay. (27) 
n {ol 
nity Expressed in terms of x, (29) can be written 
tud rays) | f1\) 1 = 
3i0ns Y 5 Ly D2 I; 5) 3] . eo (28) 
(22,)* | 2x5 2a7/) baz 
This approximates to y = x, near x, = 0 and to y 1 /6a} as v, > 0. 
It is of interest to notice that the pressure at x = 00 is given by ya} = }. 
Q A ; 
Using (6) this reduces to eae 
p 4uU R?x-* (29) 
which is the same as the expression (1) for the pressure between imper- 
on meable rollers when h 0) 
3. Distribution of pressure 
Since the function L;(€) has not been tabulated y cannot be evaluated 
9] ‘ : e . ° . by ° 
(<1 from the expression (25). For that reason the expansion of L;(€)—;(€) in 
power series was used to evaluate y. Substituting from (21) in (25) 
n ¢ e 
Ve, E cos 0 <3 )1 L»—3 ‘ 
y é e-$ sin?@ d@+-lay*. (30) 
) $7 r 
ss 0 
this 
Expanding ¢ in a power series and integrating 
(22 4 . 1) (4m-+4) 
F | ¢-€ e084 sint@ dd SY (—1)e" et =, (31) 
hea 2m! T(4m--4) 
0 m=0 = 
vers . . . 
| [he terms in the series were computed to 6 decimals as far as needed 
le! . - > . +99 
+ for é 0-2, 0-5, 1-0, 2-0, 3-0, 4-0, and 6-0. For & 6, the term in &7* 
: ; . ea 
23) | s the last one needed. An independent computation agrees within a 
| init in the last figure given 
| The computed values of y and x, are given in Table 1. 
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Fig. 3 shows the distribution of pressure. It will be seen that at 
v, = 0-29 and x, = 0°35, y is very slightly greater than 2,, so that the 
asymptotic form (27) valid for x, small is only about 3 per cent. wrong at 


“, = 0-29. The first attempt to solve (7) was made by assuming a power 
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Fic. 3. Distribution of pressure between porous rollers. 


series for y and determining the coefficients by substituting it in (7). In 
this way the following series can be derived: 


Y = x, + 623+ 300r3+ 37800713 +... (32) 


It will be seen that though the first term is the same as that derived from 
the asymptotic expansion in terms of é the series is highly divergent. 
Finding this very rapid divergence and being unable to see any physical 
reason why the solution should not be regular in the neighbourhood of 
x, = 0 an attempt was made to solve equation (7) numerically, starting 
with the values of x, and dy/dx, given by the first two terms of (32) at 
x, = 0-1, Here again the attempt failed, giving values of y which appeared 
to increase continually with increasing x,. When a descending power series 
was assumed for y the leading term was ja; and the coefficients appeared 
to be convergent, but except for values of x, greater than 1, it did not yield 
good numerical agreement with the correct values given in Table 1. 
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4, Maximum suction 
In Fig. 3 it will be seen that the maximum suction occurs at 2, 0-43 
:pproximately and that y is there 0-378 approximately. The suction is then 
T/P\3 
ts } te , pU/R\3 a 
» = (0-378 x 244)k tT Ry 0-84 (33) 
f R \k 


ind it occurs at distance 
0-43(24k R3)t = 0-95R 4) (34) 
x *$( 22H IU = ‘Jo e 
R 
from the point of contact. 


5, Straight porous band and impervious roller 

If instead of two porous rollers the problem had been concerned with a 
straight porous band running on an impervious roller (a configuration 
which occurs with the rollers used for painting walls and in paper-making 
machinery) the analysis would have yielded equation (7) after the trans- 


formation . 
3 x(96k R3)-*, y ph F 96k R3) t 
pel 


instead of (6). The maximum suction would have been 


ae eee (F 3 " 
0-378 —— (96k R*)t = 1-18—— . 3: 
eo RE) (35) 
ind it would have occurred at 
a = 0-43(96k R3)! 35k (5)! (36) 
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REFERENCES 
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THE MOTION OF A CYLINDER IN ROTATING LIQUID 
WITH GENERAL ELASTIC AND VISCOUS PROPERTIES 


By J. G. OLDROYD and R. H. THOMAS 
(Department of Mathematics, University College of Swansea) 


[Received 21 April 1955] 


SUMMARY 

A generalization is given of the theorem that any two-dimensional motion of an 
incompressible viscous liquid relative to certain axes of reference is a physically 
possible relative motion if the frame of reference is rotating uniformly with any 
constant angular velocity about an axis perpendicular to the plane of the motion, 
The theorem is shown to apply to a liquid with. quite general viscous and elastic 
properties in shear, and the way in which the force on a cylindrical obstacle varies 
with the angular velocity of the frame of reference in the general case is investigated. 
1. DEAN (1) has recently extended to the case of an incompressible viscous 
liquid some theorems first proved by Taylor (2) concerning the motion of 
a cylinder in a rotating perfect incompressible fluid. Experimental con- 
firmation of the theorems, in the case of a real liquid of small viscosity 
(water), had already been obtained by Taylor (3). It is the purpose of the 
present paper to show that the main results are of much wider applicability, 
being valid for isothermal flow of a very general type of fluid with both 
viscous and elastic properties in shear, the only restriction on the generality 
of the rheological equations of state being that the fluid is ideally incom- 
pressible and of uniform density. The essential results can be stated as 
follows: 

I. Any physically possible two-dimensional motion of an incompressibli 
fluid, defined in relation to a fixed frame of reference, is also a physically 
possible relative motion, relative to a frame of reference rotating with uniform 
angular velocity w about an axis perpendicular to the plane of the motion. 
(Rigid boundaries are to be supposed constrained to move in the same way 
relative to the frames of reference in the two cases, so that regions of flow and 
boundary conditions correspond exactly.) 

Il. The force F, on a cylindrical obstacle in a rotating tank of fluid is 
related to the force F, that it would suffer in the same relative motion if the 


tank were at rest by the equation 

F, = F,+M'[2w A R+-w A (w A R)], 
where M' is the mass of fluid displaced; here, R is the position vector of the 
centroid of the cross-section of the obstacle, referred to a frame of reference 
moving with the tank with an origin on the axis of rotation. 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 2 (1956)] 
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As an immediate consequence of II, if the cylinder is a uniform solid 
of the same density as the fluid, it will move, when unconstrained, in the 
same way relative to the rotating tank as it would move relative to a fixed 
tank, provided that the initial conditions correspond exactly. 

2. The possibility of modifying Dean’s argument so as to cover a wider 
class of materials is apparent after consideration of the restricted type of 
isothermal stress, rate-of-strain relationship which can have physical 
significance for a real liquid. Oldroyd (4) has discussed in detail the form 
which the rheological equations of state of a real material may take if they 
are to describe the intrinsic flow properties of a material element indepen- 
dently of its motion in space. The stresses in a fluid will in general depend 
on the recent strain history, so that an algebraic relation between the stress 
tensor p,, and the rate-of-strain tensor e;,, such as we have in a Newtonian 
viscous fluid, must be replaced by a set of equations involving certain types 
of differentiation and integration with respect to the time f. 

The assumption of incompressibility permits us to consider the instan- 
taneous state of stress at any point as the superposition of two stress 


systems, 1.e. , os 
Pik Pir —P %%K 


ina Cartesian system of reference), where p;,. is related to ¢;,. by a set of 
integro-differential equations of state, and p” represents an arbitrary added 
isotropic pressure which may be supposed to have no detectable effect on 
the kinematics. In general, p;, is not a deviator (ef. 4), and the variation 
of p” from place to place is determined by the equations of motion. The 
neompressibility condition implies the existence of a stream function, 
whatever the other physical properties of the fluid may be. 

We now consider a two-dimensional motion of such a general type of 
incompressible fluid to be given, the positions of particles (x,y,z), their 
velocities [u(a, y,t), v(v,y,t),0|, and accelerations | f (x, y,t). g(x, y, t), 0] 
being measured relative to certain axes Ox, Oy, Oz, without it being 
specified whether this frame of reference is fixed or moving. As Oldroyd (4) 
has pointed out, the only kinematic quantities occurring in the equations 
if state are measures of distances between neighbouring particles and the 
way these are changing at the current time ¢ and throughout previous 
history; these quantities and the material constants are not affected by 
ibsolute motion of an element in space. The equations of state therefore 
letermine the partial stress p’,(2,y,t) as the solution of a set of integro- 
differential equations whose form is quite unaffected by the fixed or 
moving nature of the frame of reference. When the initial conditions are 
specified, p,. is determined as a function of the coordinates and the time 
only. 
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Following Dean (1), we consider two motions: (1) any two-dimensional] 
motion of a general incompressible fluid, together with (2) the motion 
obtained by superposing an angular velocity w about Oz on the first motion, 
The second motion is referred to a coordinate system rotating with angular 
velocity w such that corresponding particles in the two motions have the 
same coordinates at all times. It follows, from the last paragraph, that the 
Pp; 8 are identical functions of the coordinates and the time in the two 
motions. 

Then, as in (1), the components of absolute acceleration in the two 
motions are related by 


I 


bh ~ tp) —w? 


wr, 
CX 


Ss 


B ee wy, 


> 
Jz = Jy—2w - 
. a) 
where ys, is the strearn function for the first motion. It follows immediately 
that, if the equations of motion are satisfied by one of the two motions, 
they are satisfied by the other also, provided only that 


2 ~p 7 9 9 | 9 
a 8 Pu p-—[2anp, + 4w*(a?+-y?)] 
OX Cx 
and AP2—P1) - Ps Quen, + ben?(ar?+-y?)]. 
oy Cy 


These equations are self-consistent and they determine the relationship 
between the isotropic partial pressures p” in the two motions: 


9 


Ps = pi-+2pary + bpw%(a?-+y2)+ Clb), 
where C is independent of x and y. The result I follows immediately. 

[f we now consider, in particular, flow past a cylinder in a rotating system, 
the increase in the resultant force exerted by the fluid on the cylinder due 
to the rotation can be found in exactly the same manner as for an ordinary 
viscous liquid. The proof of II is, in fact, formally identical with Dean’s, 
if p; and p3 are substituted for p, and p,, the mean pressures in the purely 
viscous case. 


3. There is nothing in the above argument which restricts the materials 
considered to be true liquids. In plastic flow of a soft solid, strains as well 
as rates of strain will occur in the equations of state, and only parts of the 
material will in general be stressed beyond the yield point. But the same 
arguments hold, provided that the material is incompressible and the yield 
criterion is unaffected by a superposed isotropic pressure. In this case some 


of the material constants may be anisotropic tensor quantities, which do 
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not occur in liquids; but then, in two motions which differ only by a super- 
nosed constant rotation, corresponding material elements are oriented in 
the same manner relative to their frames of reference at each instant, 
provided that the initial conditions are the same, and this is all that is 
required for the validity of the proofs. It is observed that even the possible 
currence of isolated cylindrical regions of plastic solid remaining in the 
elastic state and moving as rigid bodies does not affect the essential results. 
For such regions have the same density as the flowing parts of the system, 
with continuity of velocity at the yield surfaces, and it follows from II that 
no additional forces are required to maintain their relative motion un- 
changed when a uniform angular velocity is superposed on the system. 

Another extension of the theorems follows from the absence from the 
ibove arguments of any reference to the fluid considered being homo- 
geneous with regard to viscous and elastic properties; although it must be 
of uniform density, and there must be continuity of velocity at any inter- 
face between two components, in order that the proof shall be valid. In 
an inhomogeneous system of physically well-defined fluid components the 
equations of state differ from point to point, but must be the same at corre- 
sponding points in two motions which differ only by the superposition of 
uniform angular velocity; that is all that is required for the above proofs. 
The theorems I and II are therefore valid when the fluid is a mixture of 
neompressible fluids of the same density. 

One of us (R.H.T.) is glad to acknowledge the award of a grant from 


he Department of Scientific and Industrial Research. 
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THE TRANSVERSE POTENTIAL FLOW PAST A 
BODY OF REVOLUTION 


By I. J. CAMPBELL (Admiralty Research Laboratory, Teddington) 
[Received 30 June 1955] 


SUMMARY 
It is shown that in the potential flow of an incompressible inviscid fluid past a 
body of revolution set with its axis at right angles to the stream, the velocity com- 
ponents at the surface along and perpendicular to the meridians vary with azimuthal 
angle round the body in a simple manner. This is shown by entirely elementary 
considerations. 


1. WE consider the potential flow of an incompressible inviscid fluid past 
a body of revolution set with its axis at right angles to the stream (Fig. 1). 
At any point on the surface of the body the velocity may be resolved into 
the following components: 

w,, directed along the meridian. The particular value of w,, when the 

azimuthal angle @ is zero, is denoted by wf. 
W», directed perpendicular to the meridian. The particular value of w,, 

when @ is }z, is denoted by w. 
w, and w, are functions of x, the axial distance from the nose, and also of @. 

It will now be shown that w, depends on @ in an extremely simple way. 
In Fig. 2a the free stream velocity W is parallel to the meridian plane of 
zero azimuthal angle: the velocity component along the meridian at station 
x and azimuthal angle @ is proportional to W and otherwise depends only 
on x and on the angle made by the meridian plane through the position en 
the surface of the body and the meridian plane parallel to the free stream 
velocity: this velocity component may therefore be denoted by W f(z, @). In 
Fig. 2b the free stream velocity W is parallel to the meridian plane of 
azimuthal angle 26: the velocity component along the meridian at station x 
and azimuthal angle @ is Wf(a, —@) and is clearly equal to Wf(x,@). In 
Fig. 2c the free stream velocity is of magnitude 2W cos 6@ and is parallel to 
the meridian plane of azimuthal angle 6: the velocity component along the 
meridian at station « and azimuthal angle @ is 2W cos @f(x,0). Since the 
flow of Fig. 2c is simply the flow obtained by compounding the flows of 
Fig. 2a and Fig. 26 in the manner appropriate to potential flows, it follows 
that f(x,6) = f(x,0)cos@. This implies, in the notation of Fig. 1, that 
w, = wy cos 8. 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 2 (1956)] 
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The dependence of w, on @ can be derived in a similar way. When the 
free stream velocity W is parallel to the plane of zero azimuthal angle, the 
velocity component perpendicular to the meridian at station x and azi- 
nuthal angle @ may be written Wg(z, 6). When the free stream velocity W 





yW 
Fic. 1. 
Wf (x. 6) oe. 
Wf (x.6 Wf (x,-@) 2W 00s Of (x20) 
W 
W Ps 2Weos @ 
Fic. 2a. Fic. 2b. Fic. 2c. 


is parallel to the meridian plane of azimuthal angle — (7— 26) the velocity 
component perpendicular to the meridian at station x and azimuthal 
angle @ is Wg(x,7—@) and is seen to equal Wg(x, 6) by reversing the flow. 
When the free stream velocity is of magnitude 2W sin @ and is parallel to 
the meridian plane of azimuthal angle —(}7—86), the velocity component 
perpendicular to the meridian at station x and azimuthal angle @ is 


2W sin 6 g(x, 47). 
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Since the third flow is obtained by compounding the first two, it follows 


that g(x, 0) = g(x, 47)sin 8. 


This implies, in the notation of Fig. 1, that w, = w3 sin @. 

These results, that w, is proportional to cos @ and w, to sin @, are given 
in (1), where they emerge in the course of a much more detailed treatment 
of the flow past bodies of revolution than is attempted here. They seem, 
however, sufficiently important for an independent elementary proof to be 
put on record. Similar results concerning the flow over a body of revolution, 
which is rotating about a line intersecting its axis normally, emerge in (2), 
which extends the work of (1): these results also are capable of the type of 
elementary proof given here. 

[In ‘slender body theory’ (3), the flow, associated with the transverse 
component of the velocity, round any normal section of the body, is assumed 
to be the same as the two-dimensional flow round that section. This implies 
in the case of a body of revolution that w, is proportional to sin 6. 

The results are of value in picturing the flow over the forward part of a 
body of revolution when inclined at a small angle to the stream. They can 
be used to form the basis of an interpolation formula by which the pressure 
at any forward point on such a body at any (small) angle of incidence 
may be deduced from a limited number of measurements (4). They could 


also be used to deduce the streamlines over the body at any angle of 


incidence from the same limited number of pressure measurements. 
This paper is published by permission of the Admiralty. 
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SUMMARY 
ory of the subsonic motion of plane waves of finite 


is effects, without reflections and interactions. The 


applied to the problem furnishes equations suitable for 


ension of this method allows deduction of analytical 
parameters of state of a wave, due to friction; these 


numerically. An approximate method of solution of 


i. 


velocity of sound in gases ( \(yRT)). 


$4 /C). 


coethcient of friction. 


a ratio of mean to local coefficient of friction. 


absolute temperature (static). 


nensionless velocity of sound (= a/a,). 


wit.) 


dimensionless particle velocity ( V /a,). 


dimensionless pressure (= P/P,). 


dimensionless time coordinate (= agt/D). 


ince coordinate ( x/D). 
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Suffixes denote 


0 conditions in the tube before the arrival of a wave. 

01 conditions on the first elementary impulse. 

x value of parameters of state along a cross-characteristic or particle 
path. 


MW value along a characteristic (impulse path). 


1. Introduction 

THE influence of viscous friction on the propagation of plane waves of finite 
amplitude is of interest in connexion with pulsating flows. It appears that 
the recent theoretical treatments of this problem such as those due to 
Bannister and Mucklow (1), Jenny (2), and Wallace and Nassif (3) can be 
somewhat improved. Wallace and Nassif realized that a simple method 
of solution of this problem is needed for engineering purposes. In their 
analysis the assumption was made that the initial state of a wave, as well 
as all the parameters of state of a moving wave, can be connected by 
isentropic laws. Theoretical and experimental investigations made by the 
author on compression waves produced by a moving piston have shown 
that the above assumption can introduce appreciable errors. 


2. The equations of motion 

The differential equations of one-dimensional motion in pipes of a perfect, 
non-conducting gas having constant specific heats with frictional effects 
assumed proportional to dynamic pressure are 


,O(h+4V?) ah +3V?) 1eP 


J —-__~ - —_—-—_ = 0, (1) 
Cx ot p at 
a oe (2) 
Ox ot 
cP av 2 .. 
+ p—+—~ pV? = 0. 3 
ox | Ode DP i 


The above equations may be expressed in terms of variables a, V, s/R in 
the following form (4) 


9 ra y 9 a ¢ y c 
_2 aa tb Wil lé V4 O(s/R) _ (4) 
y—lévw acter ad yl ot 





2 Véa, 2 Léa, V1 a(s/R) o(s/ R) 





Bi Hon i olde Reed sos -_ (5) 
y—ladx y—lad Gx Cx ot 
9 g ’ oY a7 a 9fV2 
2 sl 41a _ a Os Ry) _ fv? (6) 
y—léx ade act y CG Da 


Equations (4), (5), and (6) represent a system of non-homogeneous quasi- 
linear partial differential equations of hyperbolic type. The mathematical 
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| lution of such a system has been studied by several authors, Courant 
ind Hilbert (5), Oswatitsch (6), and Sauer (7), inter alia. 


3, Derivation of a method of step-by-step solution 
The method of characteristics gives the following equations for the 


lirections of characteristics: 


dx 


7 V (particle path) (7) 


dx 


7 Vta (Mach wave) (8) 


dx 


7 V—a (Mach wave). (9) 


For change of parameters along characteristics the following equations are 


obtained; for the particle path 





V* 2f 4 


) d(s R) Y a he 


t; (10) 
for the Mach wave for which dx/dt = V+-a 


, 2 2fV* a 
( ac 8 ‘ se —_ > 
A OS tea oa R)+ Da aC 1) | dt; (11) 


/ / 


for the Mach wave for which dzx/dt V—a 


; 9 , 2fV3T a” 
I} = _ de hel Mh. 2 hl ~ 214 a 12 
ge ge ( ak (") 


/ 


The equations (10), (11), and (12) have been known for some time (2). 
The application of the above equations will now be explained for the 
case of a compression wave. In order to be able to determine the parameters 
of state at point 2b of Mach wave 6, as shown on Fig. 1, knowing the 
parameters of state at points 1b, P2b, and 026, a new set of equations has 
to be deduced from equations (10), (11), and (12). 
[he integration of equation (11) between points 1b and 26 gives 





2h 


2b 
| h—T, “_(a,,—a,,)+ | * d(s/R)4 ah i |e) -F| dt. 
| y J Y a 


1/ 


) 


(13) 
Ifthe points 1b and 26 are not too far apart, or if the changes in the para- 
meters of state along a characteristic dz/dt = V-+a are small compared 
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with their value, we obtain the following equation by step-by-step | and 2 


integration 






















‘i 
- . 2 ; la = : — 
Von Vip _ j (720 —Ay,) + a i(s R)s», —(s R),1»} 1 
/ } 
off V3 2b al” 
= D | \y- 1) — | j{/ap thy (14) 
m 1b 1b 
t I wet 
a Ve g° doc s V+a 
oo As CHARACTERISTIC dt 
CTERIST™ 
a er 10 a LINE/ 
= > [MAC * 
S&S . ' * +I A 
S b / | / \ar+da) , = 
\ Fé Pe / : 62k | 
& fs (06) ——¥ cross- 
cy - CHARACTERISTICS 
>, dx _,, /MACH LINE/ 
—— - : 
2a _ 
j Q 
|) on ee ns 
Saye |) eee oe ‘CHARACTERISTICS 
— dv _y /PARTICLE PATHS/ 
™ | 3 dt 
0 H T ' 
PA 
RTICLE B, B, B; oa 


The square bracket here denotes a mean value, 
equation it has been assumed that the coefficient f is constant along a 
characteristic da/dt V-+a. Equation (12) determines the change of 
parameters of state along a characteristic dx/dt = V—a (called cross- 
characteristic in case of compression waves). 
026 and 2b, Fig. 1 


Integrating it between points 
, we can write 


2h 2h 
: ; 2 fa * 2fVs a 
Von Voon y | (Hs), Agop) | y d(s R) | Da (y | ) ‘ ’ dt. 
02) 02) 
(15) 


From Fig. 1 it follows that if the particles 8, and B, are near to each other, 
then the parameters of state at point 02b and at point P2b are very nearly 
the same. The only difference between them is due to friction. On the 
other hand, the change in parameters of state between points 02b and 2b 
(dV, da,...) is appreciable because the piston has different velocities at 


the points a and 6. Hence the change of entropy between points 02) 
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y-step | and 2b is almost equal to the change of entropy of particle 8, between 
points P2b and 26. Thus, in virtue of equation (10), 


2 2b 2b 

. . r 63 2f 

d(s/R) ~ | d(s/R | y— ~ dig. (16 
| v) ; vo? D p ) 
02/ P2b P2b 


< Vd tg (Vea) dt > 
\\2b 


ra A 


Ane 





! | oy 
+ P2b 
j p<-(V+a) (dt,-dt)——> 
ATHS/ Fic. 2 
- 
‘ ‘rom Fig. 2 it follows that dt, 2dt, (17) 


where di 


denotes elementary time of propagation of a cross-characteristic. 


bias Using equations (16) and ai equation (15) becomes 





. of —" 2h 
re I ~ re lo2h " | 1) a li t ) 
J x a ((., ( . = \Y t - 2b 02b/* 
‘ j | D | J ‘tis ‘ 
> 02b 02b 
(18) 
| Adding equations (14) and (18 yields an expression for the particle velocity 
t point 24 (Fig. 1), 
E V. | , As; eo ad * 
17 a} = —4(8 R).» (s R),,} 
her, } 97 
f ty { a ~ | 
m D| | !) | j{¢2e ‘w) 
| “7 
. the ‘ 1b 
d 2l 7 
es 2h 2b 
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Subtracting equation (18) from equation (14) gives the following expression 
for the velocity of sound at point 26, 


—l,, . , Ayzt+Ag, , y—| 2h . 
Ay = : 4 (Vin—Voon) * i ye Big on [a],,1(8 R)4—(s R);»} T 
a / 


vv l j ys aad | ¢ . 
+ Ble] (e-Y—[F] \to—tw 


Y )|) « 
21 1 3 4,\ 1b 


| 
a 


2 D a J | 


02b 
The dimensionless entropy at point 26 can be calculated from equation 
(10) integrated between points P2b and 2b for a step-by-step solution, 
namely, 


2b ‘ 9 
y—1 fe rvs 2» 
y—1l =| —) 4-84 ;| \(to,—toop)- (20) 





we 


9.7 £72? 737 2b 
(8/R)o, = (8/R) po, + el (top—t poy) (21) 
D a® 

P2 

One may notice at this point that for isentropic flows equations (19) and 
(20) reduce to Vy, = Vj, and a,, = a,, respectively, as they should, because 
there is no change in the parameters of state along a characteristic 
dx/dt = V +-a in case of an isentropic compression wave. 

Equations (19), (20), and (21) enable one to calculate the parameters of 
state in the whole field of flow produced by the motion of the piston, pro- 
vided that the state at the initial (first) Mach wave and the state of 
particles at the piston face are known. 

The state at the initial Mach wave is known from the piston motion and 
an be assumed constant if the first impulse is weak. A method of calcu- 
lating the parameters of state on the face of the piston will be explained 
later in section 5. 

Once the variation of the velocity of sound and entropy along a charac- 
teristic has been found, the variation of pressure can be determined from 
equation dP 2 da 

P y—la 


When solving practical problems it is useful to introduce the following 
dimensionless quantities so that the influence of initial temperature pre- 
vailing during experiments can be eliminated: the dimensionless particle 


—d(s/R). (22) 


velocity v = V/ao, the dimensionless velocity of sound « = a/a,, and the 
dimensionless time coordinate + = a,t/D. The dimensionless distance 
coordinate then becomes x = 2/D. In the above equations a, denotes 


the velocity of sound in air in the tube before the generation of the wave. 
In the case of a rarefaction wave the method of application of equations 
(29), (11), and (12) for a step-by-step solution is similar to the one outlined 
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above, except that the impulses are now propagated along the charac- 


9 


teristics dx/dt V —a (Fig. 3). 
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Fic. 3. 
Equations for the particle velocity and velocity of sound at the point 26, 


Fig. 3, become 








, 2h 
, V, Vous Avot a ah ) » 
Me » : pe | | =~ {(s Rh)», (s R) 1») 
f (vy. 1). [2 “he si 
Dial\.\” a 
l 1 1b 
9) >} > 
loon (V2] a|” 
02h (y+1)—l= (ta, —too,), (23 
D a al ; | y | villa 
2 02b 
. Bane: Qy) Apop Y I 2b; > > ) 
2 } (Vin Voor » a 4, [a] tls Ki )op (s R)1»)} I 
Lsiy*y | a\” 
= \(Y 1)+ F (tn ty) 





\(tey—toon)- (24) 

















150 J. A. OWCZAREK 


4. Derivation of equations suitable for numerical integration for 
the change of parameters of state along the impulse paths 
For engineering applications the step-by-step method may be too 
laborious. A simple analytical method will now be derived and the general 
pattern of solutions and its trends indicated. 


(a) Compression waves 
Equation (11) determines the change in the parameters of state due to 
friction along the path of an impulse. Written in dimensionless form it 
becomes 
») 3 
Z x Vv x , 
dv yyy ; dum + ~ d(8/ R) yy + 2f ty l) | dan (25) 
np n : 


ay 
/ / 


Suffix MW indicates that changes along-a characteristic representing 
impulse path are considered; dr,,,,;- denotes elementary dimensionless time 
of propagation of impulse, Fig. 4. 

If the state of an impulse at the beginning of its motion is (v, «,...), then 
the elementary (dimensionless) time of propagation is 


dxun 

- / 92 

Ary , (26) 
vy--o& 

where dy yy dx ,;,/D represents the dimensionless distance travelled. 


From Fig. 4 it follows that d(s/R),,,, (s/R),—(s/R),. In virtue of 
equations (10) and (17) it follows that 


(s R), 4f 4 


, 
r} 


Az,+(s/R)o); (27) 


x 


Sno Nb 


where Ar, denotes the dimensionless time of propagation of the last cross- 


characteristic considered (Fig. 4). Av, can also be written as 


Ax, Ax, 
Ar, Xz Xr (23) 
V,~— a, 1.—V, 
for subsonic flow, because Ay, is negative. 
Thus, ‘ 
we al. 
(s/R)q = 4f,-y—% ——X! + (8/R)g- (29) 
r %p—Vy 
3 
tes » Ve |Axzol , ‘ 
Similarly (s/R). = 4f,y— m+ (s/R)o), (30) 
Me Lp —Vep 


where |Ay, | denotes the absolute value of the elementary dimensionless 
distance travelled by the first cross-characteristic. 


Hence, 
3 
d(s/R) sry (s/ R)a—(s/ R), “Afr 2 


Axo — |Ax, 


Xp — Vy 


(31) 
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For compression waves |Ay,, Ay,!, d(s/R)y;y- is negative; this is so 


since the paths traversed by particles decrease with distance. 
Introducing equations (31) and (26) into (25) yields 
2 Ww |\Axzo Ax, 2fv{ (y—1)—(a v) | 


dvyyy day +] : — — — _ dy yyy 
l y< N v x(v-+-a) 
(32) 


I J 


Pal (4dumw. ardanw. 7+dT7 yy ) 
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Equation (18) representing the change of the parameters of state along a 
cross-characteristic, can be written in dimensionless form as 
- or r Yr "« 
\i Aa,— 2f (y+1)+—*| Ar,. (33) 
, 1 


Combining it with equation (28) this can be re-written as 


From Fig. 4 it follows that 
A v+dvyy-—Vo1 (35) 
and Aa, = a+daygy—tg), (36) 


where the suffix 01 denotes conditions on the first impulse. As long as there 
isno shock formation, the first impulse is usually very weak and therefore 
|. Hence equations (35) and (36) reduce to 

Av, v-+dvyyy- (37) 
and Aq, r—1+day yy. (38) 
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Introducing equations (37) and (38) into (34) gives 





2 2 2f,vel(y+1)- r)] | 
dvyyt+y = - jfa—)+ 7 oy — |Axel- 
ae 1 a. x 

(39) 
For the first cross-characteristic, Fig. 4, Av, = v and Ax, = a—1, while 
Ax, = Axxo. Hence for the first cross-characteristic one can write equation 
(34) in the form 














2 -1)+(a,/v | 
y= 2 (a1) — Paral EDF (alr) ay. (40) 
; I pe et 
Subtracting equation (40) from (39) gives 
2 of, v3 (y+1)+(a,/v,) 
dvyyy = —— day + = . | (|Axzo|—|Axz|). (41) 
y—1 .( Yy Vv) 


From Fig. 4 it follows that 


(7, 7 “Te rita Te) (7,—Tq) +(ta—7p) (42) 
which for subsonic flow with vp, = 0 and ay, = 1 can be written as 
Ay lyxuw Ay 
Xxr0 ec Xun a 2 Xx! (43) 
yo, vo : ee 
From Fig. 4 it also follows that 
dxuw—|Axzo|+|Axz Xb— Xa: (44) 


Combining equations (44) and (43) gives 


IAxXz0|— |Ax.| = dxam oo (45) 


f 

LI 
1+{1/(a,—v,)} 
As an approximation, the mean values of the parameters of state along 
cross-characteristics and particle paths can be assumed to be equal to the 
mean values of the parameters between the first and last impulse con- 
sidered. Thus 
sai . (46) 
r 3 


and V, : (47) 


— 


aS Vp, = 0 and ag, l. 


From equations (41), (45), (46), and (47) it follows that 


2 , Pl (y+1)jyt+at1](v-+a—1) 
dv =: a = vy: 6 OEE ia eee. Nee a SE " (48) 
at y—1 uwtta (a-+-1)(a—v-+-3)(v-+a) XMW 
Introducing equations (46) and (47) into (32) gives 
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Dropping suffix MW and denoting the ratio of the meanj to local coefficient 


f friction f,,/f by k, the following equations for the frictional changes along 


the path of an impulse (characteristic dx/dt = V-+-a) can be written 











2 v7 (y+1)jyptatl1](v+a—1) , 
dv da Lf’ ss, 1 50 
= y—1 itl os 
ind 
2 2v* ((y—1)v—a 2va(vta—1) | " 
ly da , nh, Na a 
y—1 l x | x (a v-+3)(a-+1)24 °* _— 


Equations (50) and (51) can be integrated numerically. In general, the 


nitial values of the parameters of state v and a cannot be connected by 
sentropic laws, but have to be determined from equatiens given in section 5. 
Elimination of parameter y from equations (50) and (51) yields the following 
lifferential equation combining v and «a: 
f__ Ff x — (y—1)v ]|(a—v+3)(a+1)?— 
2 theva?(s r—1)+-ka(a- L)[(y-4 l)v- x+1](v+-a—1)} 


ka(a- lf (y- ljyvta+1](v-+a—1)} 


(52) 

In the subsonic region above the line v 1—a, which is the main region 

f physical application dv/da < 2/(y—1); dv/da = 2/(y—1) for a = 0 and 
] x. In the region below the line v l—a, dv/da > 2 (y- 1). 


Fig. 5 shows the pattern of solution of equation (52). 
Combining equations (50) and (51) yields the following relationships 


theva(s r— 1) kl (y + 1)yvta- 1 |(v 1 qy— ] 





} dy, (53) 


(a+1)*(a—v+3) (a+1)(a—v+3) 
(“| ; 1) y= | Z| a—(y Lv] 
2 2 ( x) | x 
to me 1) kl (y | 1 )v 1 oy +1 ](v , le, (54) 
(a+1)*(a—v+3) (a+1)(a—v+3) | 
In the subsonic region di 0 and dx < 0. Hence the effect of friction is 


to reduce both the particle velocity and the temperature. This tendency 
s indicated on Fig. 5 by arrows. 

Expression for the frictional pressure change along the path of an impulse 
s obtained from equations (11) and (12) after the velocity of sound and 
ntropy terms have been eliminated, and the pressure term introduced with 


Mean along a cross-characteristic or a particle path. 











154 J. A. 


OWCZAREK 
the aid of equation (22). An analysis similar to the one outlined on the 
last few pages gives the following equation 
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where 7 = P/P, = P/P), represents the initial dimensionless pressure of 


an impulse. In the derivation of equation (55) it has been assumed that 
1), by analogy with equations (46) and (47). 


a L(x ! 
TT 5 (71 


It may also be 
noted that equation (17) has not been used in the derivation of equation 
(55). Equation (17) is only an approximation when considering points far 
apart on a cross-characteristic. 

As a first approximation, solutions of equations (53), (54), (55) are arrived 
at by assuming constant values of functions on the right-hand side. In most 
engineering applications, the changes in the parameters of state are small 


enough to justify the above assumption. 


(b) Rarefaction waves 

The method of derivation of equations corresponding to (53), (54), and 
(55) in this case is similar to the one described above. The main difference 
lies in the fact that the rarefaction impulses propagate along characteristics 
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Equation (12) becomes 


rs x vs x wis 
dvyry Axyyy5 d(s/ R)yny—2f (y—1)+ dt yy: (56) 
y l Y t Vv 
[he dimensionless time of motion of an impulse is 
dxww ~— 
At yy ° (97) 
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From Fig. 6 it follows that d(s R) am (s R), (s R)... By analogy with 
what was said for compression waves 
3 
. vw Ay,—Ay, 7 
d(s R un 4fy = Xa é Xr0 (58) 
My Vp Oy 
Here d(s RB) ary O as Ax Axo 
Equation (56) combined with equations (57) and (58) takes form 
2 Xl 3 Ay ' Ax “() v x \dyww 
dv ary dayy—4f,— = = 2f—|(y—1)+- —, 
I v.20 Vp +a, x v| a—v 
(59) 


When written in dimensionless form, in which the entropy term has been 
modified as in equation (18) for compression waves, equation (11) for cross- 


characteristic da/dt V +a takes the form 


z..., (60) 
x 


Ax 
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Using equations (37) and (38) we obtain for the last cross-characteristic 
(Fig. 6), 








9 9 3 A 
~ “ Lor Ya me; Xx 
v+dyyy = — 7 (a—1)— day + 2f, J+) | os 
: a : a Xz Vea hy 
(61) 
For the first cross-characteristic equation (60) becomes 
9 3 
2 wv x] Ay.. 
y= ———— (a—1) 4 2f, 21 (y+ 1) — 2 | — (62) 
o Xz Vy Ya Xy 
Subtracting equation (62) from (61) yields 
9 3 1A is nt . 
«“ of Ya Xe | OX2 Xx0 +6 
dvygy = ——— day + 2f2—| (y+1) —= |-_*—_—= . (63) 
ia l Xp v _ Vy Xz 
From Fig. 6 it follows that 
i Ta) om (ta—7-) (7, —7,)-+ (74- Tp) (64) 
which can be written as 
Ax dy ww A 
x0 | Xuu | Xe +» 
+: X»— Xal +X. (63 
Vo OX, Y—vy Vn-{- Oe 
From Fig. 6 it also follows that 
Axet ldxaw|—Ax20 Xb— Xa'- (66) 


From equations (65) and (66) we obtain the relationship between 
Ay,—Ay,9 and |dy4,;,|.. This relationship and equations (46) and (47) intro- 
duced into equations (59) and (63) with suffix MW dropped and the ratio 
f,/f denoted by k, yield the following equations 


dv . da f = (y= ress as - 7 | dx’, (67) 





y—1 “a—v{ x © (vta+t3)(a+1)? 
a. a rp CLs 1)v—a—1](1+v—a) dy. (68) 
y—1 ; (a+ 1)(v+-a+3)(a—v) ‘ 


In the above equations dv and da represent elementary changes of para 





meters v and « due to friction along the path of a rarefaction impulse, and 
dy| represents the absolute value of elementary dimensionless distance 
travelled by it. Equations (67) and (68) can be integrated numerically. 
The initial values of parameters of state (v,«,...) cannot, in general, be 
connected by isentropic laws. 


Elimination of parameter y from equations (67) and (68) gives the 


following relationship connecting v and « 


{—2[ (y—1)v+a](a+1)?(vt-a+3)—4hva?(1+-v—a)- 





dv 2 ka(a- L)fa t+1—(y- lv] 1 y—«a)} 
dx y—1 {—2(y—1 pv faf(at+1)2(v+a+3)—4hvaX(1+-r—a)+ 
t-ka(at1)fa+1—(y+1)v](1+rv—a)} 

(69) 
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Fig. 7 shows the pattern of solution of equation (69). Along lines 
v—1 and (y+1)v = a+1, dv/da —2/(y—1). In the regions 1 and 4, 
Vv=a-1 (y+1) v-a-1=0 
x ? p 
i 
V4 
Py 
¢/ 
Ww 
1-0 —T 
} | 
ok 
: —_—_ 
0 , 1 0 v 
Fic. 7 
dy/dx - 2/(y—1); in the regions 2 and 3, dv/da > 2/(y—1). From 


? equations (67) and (68) the following relationships follow : 


di i I 2| r+(y 1 )v| ; 
2 a—v\ \ 
_ _Akva(i+v—a) _ Ma+l—(y+l\(l+v—a)) dx|, (70) 
(a+ 1)?(v+-a+3) (a+1)(v-+a+3) | 
rm i 4 1 v® (2la+() -1)v] ” 


2 2 XY v\ ¥ 
thva( 1+ x) ki r+ i —(y- Dale! +-y—a)| 
(a 1)*(; r+3) (a+ 1)(v-+a+3) 


> Inthe region of applicability (subsonic flow), dx > 0, dv < 0. Hence the 
effect of friction on rarefaction waves is to increase the temperatu.e and 
decrease the particle velocity. This tendency is indicated on Fig. 7 by 
arrows, 

An analysis similar to the one outlined above was performed using 
equations (11) and (12) with the pressure term in place of the velocity of 
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sound and entropy terms. It gives the following expression for the frictiona| 
pressure change along the path of a rarefaction impulse 


pany Se ee 
" (24+1/a)7+1 «a a—v| x 
k| r+ 1—(y l)v|(1 1-y—c)) 
- (a+-1)(v-+a+3) | 
Equation (72) indicates that friction tends to increase the pressure of a 
moving rarefaction wave. For v—> a (M > 1), both y and (a—v) tend to 
zero. Hence, when performing calculations near M = 1, the dimensionless 
time coordinate should replace the distance coordinate given by equation 
(57). 
Similarly as in case of compression waves, approximate solutions of 
equations (70), (71), and (72) are obtained by assuming the right-hand 
side functions to be constant. 


5. Methods of calculation of the initial state of a wave 
(a) Compression waves 

In case of compression waves generated by a moving piston of known 
velocity as a function of ¢ (Fig. 8), the velocity of sound can be calculated 
from equation (18). 

Using the notation of Fig. 8, equation (18) can be written in the follow- 
ing dimensionless form 

l y3 


x X14 a (v—V9}) Ply 1) “tor ee 


“~ x 


(?—Tg,)- ( 





As a first approximation for the calculations, «, should be taken as a mean 
value between a», and a;., «;, denoting the isentropic value of the velocity 
of sound corresponding to the piston velocity v. For weak initial impulse, 
one can assume vp, = Oand a, = 1. The specific entropy can be calculated 
from equation (10) by integrating between 7 = 0 and 7. The pressure can 
be calculated from equation (22), once « and s/R are known. 

In case of waves generated by efflux of gas from a vessel into a tube, we 
may expect that the assumption that point (v, «,z,...), Fig. 8, is at the tube 
inlet will be justified. Then, if 7 and « are known, the particle velocity can 
be determined from the equation for a cross-characteristic, namely, 
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(b) Rarefaction Waves 
In cases of waves generated by a moving piston, Fig. 3, the particle 


velocity is known from the piston motion. The parameters of state which 
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OF 





have to be calculated are the velocity of sound and the specific entropy. 


The pressure can be determined from equation (22) once a and s/R are 









known. 
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The velocity of sound can be calculated from equation (60) written in 
he for) 
Xo] > J Vo) TAY 1) (y 1) : (rt T 91): (79) 
2 x, , 





[he specific entropy can be calculated from equation (10) by integrating 





petween 7 {) and a and IS 
ps 
) , . » 7 - 
R (s/R)o, +2f,y—r. (76) 
VP 
In case of waves generated by a sudden discharge of gas from a tube 
Fig. 9), the pressure at point p is equal to the external pressure. The 
particle velocity can be calculated from the equation for a cross-charac- 
teristic dax/dt V-+-a, which when combined with equation (22) takes the 
rm 
"ee a py 
01 1 
; . 2). l(y | v.— a, |(T—T9;)- (77) 
jj x, 
The value of the specific entropy can be calculated from equation (10) inte- 
grated between 7,, and rt. The value of the velocity of sound can then be 


determined trom equation (22). 


In case of waves generated by discharge of gas from a tube through a 
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vadve (Fig. 10), the values of the parameters of state at point p should 
satisfy (a) the equation for the particle path (10), (6) the equation for the 
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cross-characteristic (11), (c) the equation for expansion to the external 
pressure corresponding to the exit area determined by the position of the 
exhaust valve. Isentropic position diagram should be used for the purpose 
of integration. We may expect that the assumption of isentropic steady- 
flow expansion through the short valve opening will be sufficiently accurate 
for engineering purposes. 
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Concluding summary 
The method of characteristics furnishes 
the < 
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the following equations for the 


iccoustic velocity, and the dimensionless entropy for 
For compression waves we have 
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pressure. For compression waves we have 
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waves are different from those corresponding to the rarefaction waves. 


From the finite difference treatment the following differential equations 


are obtained for the particle velocity, the accoustic velocity, and the 


(54) 


(55) 


(70) 


(71) 


(72) 


In the deduction of the above equations it has been assumed that the 
particle velocity and the acoustic velocity along a cross-characteristic are 
constant. Their value has been taken as an arithmetic mean between thie 
initial and the final states 

In the case of a pressure wave having a peak, the integration should be 
performed in two steps, the first terminating at the intercept with the 
path of the wave peak, and the second extending from this intercept to 


the required point. This is so because equations describing compression 
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From equations (53), (54), (55) and (70), (71), 72), it follows that the effect 
of friction on a compression wave is to reduce the particle velocity, the 
temperature, and the pressure, while the frictional effect on a rarefaction 
wave results in an increased temperature and pressure and a decreased 


particle velocity. 
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MEASUREMENTS OF ATTENUATION OF COMPRESSION 
WAVES OF FINITE AMPLITUDE IN AIR AND 
EVALUATION OF THE COEFFICIENT OF FRICTION 


By J. A. OWCZAREK (Queen Mary College, London) 
[Received 10 September 1954; revise received 14 April 1955] 


SUMMARY 


The purpose of this paper is the investigation of the influence of viscous friction 
on plane non-steep compression waves generated in air in a tube by means of a 
pneumatically operated piston. 

The tests described include four series of experiments. In every series the shape 
of the pressure wave (pressure—time curve) was determined by means of a cathode-ray 
oscillograph using ‘ variable capacity’ pressure gauges in four test stations 10 ft. apart. 
Also the displacement—time curves of the piston generating the waves were recorded 
for the purpose of determining the piston velocity curves. Each series of tests has 
a different piston velocity curve. 

The acceleration of the piston varied over the range 7,260 to 27,200 ft./sec.?, and 
the highest pressure of the generated waves varied from 6-5 to 8-7 lb./in.? gauge. 

The measured pressure waves show a decrease in pressure of the generated waves 
with distance travelled and an increase in the time of duration when compared with 
isentropic waves. 

The equations for calculating the coefficient of friction in non-steady compressive 
motion of a fluid are presented and the calculated and correlated values of the 
coefficient of friction plotted. 

It has been shown that the coefficient of friction depends on the state of the 
boundary layer between the sections considered and that it must, therefore, depend 
on the Reynolds number, Re,, based on the distance travelled by the mean particle 
between these sections. The mean particle, located somewhere in the middle between 
the sections considered, should be a fair representative of all the particles of fluid in 
motion between these sections. Thus the assumption prevailing that steady flow 
coefficients corresponding to Reynolds number based on diameter should be used 
in the case of non-steady flow has been disproved. A comparison has been made 
between the coefficients for steady flow (with developed boundary layer), and those 


for non-steady flow. It has been shown that the acceleration in’ eases the values of 


the coefficient of friction. 


Notation 
A_ cross-sectional area or acceleration. 
a velocity of sound (= v(yR7)). 
a, velocity of sound before generation of pressure wave. 
(perimeter. 
D hydraulic diameter (= 44/C). 
f local coefficient of friction (= 27/pV*). 


(Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 2 (1956)] 
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L characteristic length l ft. 

P pressure. 

P, pressure before generation of pressure wave. 
R gas constant. 

gs specific entropy. 

t time. 

T absolute temperature (static). 

V particle velocity. 

x distance. 

x dimensionless velocity of sound (= a/a,). 
ratio of specific heats (= c,/c,). 

[i absolute viscosity. 

dimensionless velocity (= V /ag). 

7 dimensionless pressure (= P/P,). 

po density. 

o Langhaar’s parameter defined by equation (24). 
7 shear stress or dimensionless time (= a,t/L). 


y dimensionless distance (= 2/L). 


Suffixes denote: 
0) initial conditions. 
01 conditions on the first impulse. 
r mean value along a cross-characteristic or along a particle path. 


f frictional influence along a characteristic. 


1. Description of apparatus 

THE apparatus shown in Fig. 1 was originally built by J. S. Glass.7 It 
consisted of a steel tube, 2 in. internal diameter, 51 ft. long, a pressure 
vessel, a piston, a piston releasing mechanism, an electronic indicator, and 
acompressor. The tube, into which the piston was propelled, was polished 
inside along its whole length and had four pressure measurement stations 
rbout 10 ft. apart. The first test station was 98 in. from the piston releasing 
mechanism. The last station was 145 in. from the open end of the tube, so 
that a measurement of the whole pressure wave at this station could be made 
before the arrival of the reflected wave. At the end of the tube there was 
in orifice }} in. diameter. Its purpose was to let part of the generated wave 
escape from the tube so that the strength of the reflected wave was just 
enough to stop the piston half-way along the tube. At the end of the tube, 
near the releasing mechanism, there were twelve electric contacts. The 


piston was 48-64 in. long and was composed of two aluminium 2-in, diameter 
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pistons with rings, connected by an aluminium alloy tube 1} in. diameter, 
It had a flat wave-producing face thus ensuring production of a plane wave, 
Before release, the rings at the front of the piston touched the 12th contact. 
while the rear end of the piston was sealed off from the releasing mechanism 
and the pressure vessel by means of a rubber membrane. The purpose of 
the slots in the tube was to let the compressed air propelling the piston and 


Test station IV Test station 


ll Test station Il Test station 


——— 145-1" 129-1"- 1203" 


a. ae —- . 





z 4 5 6 7 Q 9 ry 
Fic. 1. Layout of apparatus. 

1. Orifice. 6. Air compressor.” 

2. Stand. 7. Electric contact. 

3. Condenser pressure gauge seating. 8. Piston releasing mechanism. 

4. Tube, 9. Pressure gauge 

5. Muff. 10. Pressure vessel. 


leaking through the rings of the propelled end (rear) of the piston while it 
was in motion escape into the atmosphere without interfering with the 
wave producing end (front) of the piston. As the piston started its motion, 
the rings of the wave producing end stopped touching the 12th contact. 
This event was at once recorded on the time marking signals of the indicator. 
The rings of the propelled end of the piston would now start touching all 
the twelve contacts on its way and send signals to the time marker of the 
indicator. In this way the time taken for the piston to travel through the 
distance equal to its length was recorded. From the known distances 
between the contacts the velocity of the piston could be determined as a 
function of distance. At the propelling end of the piston there was a catch 
through which the piston was held in position by the plunger of the releasing 
mechanism. A channel in the steel body of the releasing mechanism con- 
necting the high pressure vessel with the piston of the plunger had a plug 
which obstructed the channel in the normal position. When the plug, 
operated by an electromagnet, was pushed out of its normal position, the 


air from the high pressure vessel acted through the now opened channel on 


the piston in the releasing mechanism, moving it downward and thus 
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unhooking the main piston in the tube. Now the compressed air, acting 
on the propelling side of the piston in the tube would break the rubber 
membrane and accelerate the piston. The transient pressure recording 
equipment was a single channel electronic indicator utilizing condenser 


pressure gauges and a rotating drum camera. 


2. Experimental results 

The experiments consisted of four series of tests on different compression 
waves measured at four test stations along the tube. About seven photo- 
sraphs of pressure waves were taken in each series at each test station. 
Figs. 2, 3, 4, and 5 show the experimental and calculated isentropic pressure 


waves for series 2 plotted in dimensionless coordinates 
7— | (P—P,)/P, = AP/Py and +r=a,t/L. 


The dimensionless time coordinate 7 denotes the time which has elapsed 
from the beginning of the piston motion. Comparison of Figs. 2—5 indicates 
that the pressure decreases substantially with distance and that the time of 
duration of the waves increases with distance as compared with the time 
of duration of isentropic waves. This effect is indicated by the horizontal 
extension of the plotted values beyond the point 12 on the isentropic 
line 

The maximum pressure ratios (attained pressure to atmospheric pressure) 
nseries 1 2, 3, and 4 were 1-444, 1-478, 1-546, and 1-592 respectively. The 
maximum velocities were 304, 324, 2376, and 401 ft./see. The maximum 


ccelerations in ft./sec.2 were 16,800, 19,000, 26,100, and 27,200. On Figs. 
2 and 3 points calculated using a step-by-step method (1) have been plotted. 
In the calculations, coefficients of friction equal to those corresponding to 
steady flow with a developed boundary layer for the given Reynolds 
number, R pVD pu, have been used. The calculated pressure wave, 
corresponding to the 


given curve of the piston velocity, lies above the 
isentropic wave whereas the experimental results lie below it. This dis- 
crepancy in the positions of the curves can be attributed to the effect of 
heat transfer to the tube and to the face of the wave producing piston, 
ind to the non-one-dimensional flow conditions. 

If the calculated values for station I are moved so that they correspond 
to the experimental curve (‘corrected values’ of Fig. 2), and the calculated 
points for station II are shifted by the same amount, then the calculated 
values are below the experimental ones at this station. This effect is 
especially pronounced in the initial part of the wave, thus indicating that 
the coefficient of friction in the initial part of the wave has been taken too 


high and in the final part too low. 
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In (1) the effect of heat transfer on the experimental results has been 
nvestigated. The conclusion reached was that for a moving wave of 
pressure ratio up to 1-6 produced by a piston in a tube which has a tem- 


perature equal to that of air before compression, the influence of friction 


s predominant. Hence difference in the change in the shape of the experi- 

mental waves and the calculated isentropic ones could be attributed to the 
effect of friction only. For the shape of the waves produced on the face 
f the piston the effects of friction and heat transfer were found to be of 


the same order of magnitude. 


3. Derivation of formulae and evaluation of the coefficients of 
friction 
The equation determining the change in the parameters of state due to 


friction along the path of a compressive impulse can be written in dimension- 
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less form with finite differences as, (1, 2), 

Ay, = IAG OF eft] ae a) 

™MW> 
1 —"y [P] DIL a 

because the changes in the parameters a, v, and P are small. Suffix f has 
been introduced into equation (1) in order to distinguish frictional from 
other changes. A square bracket denotes a mean value along a charac- 
teristic, Ary; = A Aty,y-/L represents dimensionless time of propagation 
of the impulse, and can be written as 


Axa 
At iny “. (2) 
[p)+[a] 
Introducing equation (2) into (1) gives 
x| AP, 2 —(y—1) v/a}! 
aie "I PI- D f [Y’ yu ied Axanr- (3) 
y [P] 4 [»]+[e] 
The equation determining the change of the parameters of state along a 
cross-characteristic, written in dimensionless form with finite differences, 
can be written as (1, 2) 
-AP, Ff, 9 {1+(y—1)(v,/a,)} 


> ,—v 


I a si 
where Ax, (v.—a,)Ar,. 
Suffix 2 denotes a mean value of parameters of state between the first and 
last characteristic considered. Ay, is the dimensionless distance covered 
by a cross-characteristic. 
From Fig. 6 it follows that for the last cross-characteristic 


Av. vo Av, Vo1 (5) 


and AP, = P+AP,—Py, . (6) 


where v and P are the dimensionless particle velocity and pressure at the 
beginning of motion of the impulse considered, Av, and AP; are the fric- 
tional increments, vj), and P,, are the parameters of the first impulse 
produced and are assumed constant. 


For the first cross-characteristic, Fig. 6, 


Av, V—Vo}) (7) 
AP. P—P,, (3) 
and Ax, Ax,0- (9) 


Hence for the first cross-characteristic equation (4) can be written as 


,P—fy 2f, etit(y—1)(v,/«,)} 


v—v . ve 
01 y FP. D/L” Xe—vV 


J Ir 


Axzo- (10) 
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For the last cross-characteristic equation (4) becomes 





Xy AP, Xy P~ie , 2f, {1 T (y—1)(v,/a2)5 Ax 





, Av,z+(v—vp)) + a : 
f l 4 zx 
P  e » D/L * p—Vz 
(11) 
Subtraction of equation (10) from (11) gives 
xr AP, fp 2tltly—VW(r2/%2)5 9 
Av; » P Dil Vy eee ape  (Ax,—Axz0): (12) 
N / I a r x 
+ 
- ayy, _ @rAvpersheyy) 
_Aot| As TR 
T\2 pone a CROSS 
1\g CTE , 4 \ CHARACTERISTIC 
¥ THARAY an \ 
(v,@) Gz | 
Rs cross Se | (Voss ®o1) 
x ‘ | ow 01 
> x ~ Sit 
se 
&/% 
SY fe 
9) 
= 
g 
« 











In compression waves without reflection or interaction Az, is positive and 
| Ay, negative, Ay,» has a larger negative value than Ay,. Hence, with the 
absolute value of the difference between Ay,,. and Ay, equation (12) becomes 
x, AP, 2]. (1-4 (y—1)(v, Xx) 5 A 


9 
Vy 


| (ro AXz|: 13 
P. D/L 1,—Vy Xxr0 Xx (13) 


Av; 
Elimination of Av, from equation (13) using equation (3) gives the 
| following equation 
AP, [a]+a,([P\/P,) D 
) 


P| Dy L na 





—kv, - : =F : Axz0 Ax, 


p)2=o—Db/aD Val 4 


L+[r}/Lo] 
in which the ratio of the mean value of the coefficient of friction along a 
cross-characteristic to the value of the coefficient along the characteristic 


considered is denoted by f,/f = k. 
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In addition to equation (14) another equation for the determination of 
the coefficient of friction from the experimental results can be deduced, 









As Aty = Atyyy—Artyyyig (Fig. 7), where Ary))-;, denotes the time taken 
~Gef {_! (v+Avy,a+Aar) 
3 1 og Bll f 
<1 | 
: ADIABATIC AG, 0) 
bt | 
| | 
| | dy | 
MW —=V+t+a 
dr 
ISENTROPIC 








|; ——- > 
AX uw x= 
L 
Fic. 7. 
for an isentropic impulse to cover distance Ay4,,,, equation (2) can be re- 
written in the form 
Axwn 4 Av; | ae (15) 
_ > ») » 
Ar,+ Artanis 2 2 
as [vy] = v4 (Av,/2) and [a] x+(Aa,/2). The equation combining the 
velocity of sound, pressure, and entropy can be written in dimensionless 
form as 
Aw y—1, ,AP.. y—1 
f= ~_ [a] —f4 "  A(e/R). (16) 
2 fy [PI] 4y 
Introducing equations (16) and (3) into (15) gives 


3—y, ,AF, Ui 2111 —(y—1)[v/a]} 
dy lol Pp) D AL | [v]+[a] XMu 


va Bist, - 
} A(s/R) Xu . (17) 
4y Ar; At awiris 
The change of entropy along a characteristic can be represented by the 
following equation (1) 
of p 
A(s R) 1 ng 2 Aa (18) 
D/L’ oo” ”” 
rt 
where Az,, represents the difference in the time of movement of the particles 
at the limits of the distance covered by the impulse. Suffix a2 denotes a 
mean value of the parameters between the first and last characteristic 
considered. 
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tion of Introducing equation (18) into (17) and rearranging gives 
gquced 
luced. (1 ; | v1) 3 F AP, 
ti 9 
iken »  L+(Atapis/AT,) 1, P|D. y-—!1 f ve [v] | «| v, AT, 
ro] § 7 T 9 42 27,2 a re a 
1 l—(y—1)| v/a \ L 2 x2 [v2 |{1—(y—1)[v/a]} Anaya 
l SAM 
¥ l x | 
(19) 


The second term of equation (19) is negligibly small. For example, the 
/ error in the value of f introduced by neglecting this term im calculations 
for series 2 between test stations I and IV was found to be about 0-2 per 
ent. Hence, to a good approximation, 

v)(1/{«]) 3—y AP, 
AT Wiis Ar) fy |P| D 
) y*| ,1—(y l)|v LL 
[y]+[o] 


From experiments the frictional pressure drop AP,/P, was known for given 





(20) 





x 
ro 


\yy- and also the frictional time increase Az, as a function of P/P, for a 





‘ope. | given Ayyyy- Hence the coefficient of friction f could be determined from 
, equation (20) for each series of tests. As Av, could not be determined with 
great accuracy, it is expected that the coefficients of friction calculated 
5 ising equation (14) are more reliable. Using equation (14) the values of 
the coefficients of friction have been calculated and compared with those 
- the obtained using equation (20). In the calculations the mean particle velocity 
nless for a cross-characteristic has been determined from the equation 
@ l }((v] Vo1)> 
the suffix 01 denoting conditions on the first characteristic; the mean 
velocity of sound has been determined from the equation a, = }({«]+-a 9). 
| The value of « has been taken as the isentropic value corresponding to the 
| given piston velocity, because the cooling effect of the piston and the tube 
| walls reduces the value of « almost to the isentropic value. For example 
7 for characteristic 7 in series 2, « 1-0445 while a. 1-0450 and the value 
fx calculated with frictional effects with the ‘steady flow’ coefficient of 
the 


lriction was 1-0466, The pressure values [P| and P, = }({P]+P,) have 
been taken from the experiments. The coefficients of friction calculated 
1s) | forall four series of tests are plotted separately on Figs. 8, 9, 10, and 11, 


ind together on Fig. 12. The values of the coefficients have been calculated 


les for experimental results obtained between test stations I and II, I and 
sa (II, and I and IV. In Fig. 12 the values of the coefficients of friction for 
sti laminar and turbulent boundary layers in steady flow along a flat plate 


have been given for comparison. 
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When considering the experimental results the distance x in the Reynolds 
number denotes the distance travelled by the mean particle between the 
test stations considered. In the flat plate equations the x coordinate 
denotes the distance from the leading edge of the plate. The values of 
the coefficient of friction could not have been calculated from the very 
beginning of the pressure waves because the beginning of the shock forma- 
tion affected this region. This shock formation effect is especially pro- 
nounced in series 4. Thus on Figs. 11 and 12 the first calculated values 
of the coefficient of friction for this series correspond to the region where 
the transition boundary layer has already been formed. The values of the 
coefficient of friction shown in Fig. 11 in a dotted rectangle appear to be 


already affected by the shock formation and are therefore too low. The 


curves in Figs. 8-11 show a decrease in the value of the coefficient ol 


friction when the experimental results used are those between test stations 
I and II, I and III, and I and IV. As the acceleration and the pressure 
gradient for a given impulse path do not change appreciably, the explana- 
tion is that the value of the coefficient of friction depends on a parameter 


which incorporates the distance travelled by a particle. For correlation 01 


the coefficients of friction in the laminar boundary layer region the correla- 
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tion equation as derived by Moore (3) for a flat plate has been used. It can 


be written in the form 





rA x2 dA 
(( 14- 2-555 — — 1-414: - | 21) 
} F , V2 V3 dt ' } ; 
——n 
f 
seal ‘. , ae (MOORE) J tre = — 
UL 47.c0c XS we at 
i 3 73 414 : 1+712 (24) 


ee LAMINAR B.L TRANSITION B.L 
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% St - 
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In the above equation f., denotes the coefficient of friction for steady flow, 
the coefficient in an accelerated flow, and A dV /dt the acceleration. 


lenotes the distance from the leading edge of a plate; in our calcu- 
tions it becomes the distance travelled by the mean particle between the 

test stations considered. The correlated coefficient of friction for the 
minar boundary layer is shown in Fig. 13 and can be represented by the 

lat n ~ 

wo 0-025 


fin ae 22) 


(Re .) a6 


lt is somewhat above the line corresponding to the laminar boundary layer 
fora flat plate in steady flow. It is, however, an established fact that the 

efficients of friction in tubes, calculated for mean fluid properties at a 
given cross-section, are larger than those for flat plates which are calculated 
lor the undisturbed fluid properties. One of the first to notice this fact was 
Wieghardt (4 





Eckert (5) has shown that the results of investigations on 
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the friction in pipes can be applied to flat plate problems if in the resistance } tion 1 
formula the values of the parameters of state corresponding to the centre | The c 
of the tube are used. He has found that the ratio of the coefficient of 
friction for fully developed turbulent boundary layer using the centre values 
to the coefficient of friction using mean values is 0-7. Apart from steady 
flow along a flat plate, another comparable limiting case of non-steady flow | The ¢ 
in pipes is the steady flow in the ‘inlet length’ of a pipe. This type of flow | (26) is 
has been investigated by Langhaar (6). In an attempt to make use of 
his findings the following definition of the coefficient of friction will be 
assumed 


' D -df It is ¢ 
I aes (23) 

4 toV2dx to tu 

frictic 


As shown by Elser (7), the coefficient of friction defined by equation 


(23) depends on the Mach number as well as the Reynolds number. For _— 
small Mach numbers, however, this effect can be neglected. Introducing a 
Langhaar’s coordinate — 
- 4xr/D 4(a/D)? (24) betwe 

Re, Re, 7 Serie 

and noting that dP/dx = —d(AP,)/dx and that in the ‘inlet length’ in 4 
steady flow the changes in the dynamic pressure and the Reynolds number, . 
Re,, are small, one can write Serve 
p_ 1 MAP| SeV2) is . 

‘ Re, da , 

The value of AP,/}pV? is presented by Langhaar as a function of o — 
in his Table 2. As the values of o for the experiments performed are : 


rather small, a slightly extrapolated value of AP,/}pV* had to be used. = 
The result which has been obtained by applying equation (25) with Re, | °@™ 
representing a mean value of the Reynolds number between the charac- Te 
teristics considered, is shown on Fig. 13. It indicates the usefulness of t 
Moore’s correlation equation for non-steady flows in tubes with laminar [| Fror 
boundary layer. tion 

From inspection of Fig. 11, which shows the distribution of the coefficient | _ pres 
of friction for series 4, it is seen that the turbulent boundary layer | byn 


has almost been reached. Bearing this in mind, and realizing that the | x,,,, 
correlated values of the coefficient of friction should be somewhat above } Rey 
those for the turbulent boundary layer on a flat plate in steady flow, a It 


new correlation equation has been introduced for the transition region. | serie 
Moore’s correlation equation for laminar boundary layer applied to the | dev 
region of transition makes the values of the correlated coefficient too low. | initi 
There is no theoretical justification for the correlation used for the transi- | com 
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tion region; there is, however, no known theory applicable to that region. 


[he correlation equation used for the region of transition is 


At -712(F5) |. (26) 


The correlated coefficient of friction in the transition region using equation 
26) is shown on Fig. 13, and can be approximated by the equation 


0-092 14500 


(Re,)* - Re. 


I 


ft frre (27) 
It is difficult to draw many conclusions about transition from the laminar 
to turbulent boundary layer from the distribution of the coefficient of 
friction because of the one-dimensional analysis employed. Several para- 
meters based on mean fluid properties and corresponding to the transition 
points determined from Figs. 8-11 have been calculated and are given 
below. The calculation of parameters has been done for mean particles 


between the indicated test stations for one point of each series of tests. 


‘é rie 8 ] 


Test stations I-II. 2,., 2-0 ft. 

Re 3°47 x 106; M V/a = 0-198; P/P, = 1-323; xA/V? = 0-447, 
Series 2 

Test stations I-IV. 2,,.., 1-82 &&. 


Re 3°56 x 108; M V /a = 0-218; P/P, 1-362; zA/V? = 0-411. 


est stations I-III. a%,.,, 1-66 ft. 

Re, = 3°31 10°; M V /a = 0-224; P/P, = 1-351; 2A/V? = 0-486. 
Series 4 

Test stations I-IV. 2;,., 1-13 ft. 

Re 2:4 10°; M V /a = 0-236; P/P, = 1-366; xA/V2 = 0-384. 


From the values of the parameters given above one notes that the transi- 
tion of the boundary layer appears to begin in all series of tests when a 
pressure ratio of about 1-35 has been reached, while the distances travelled 
by mean particles from the beginning of the motion to the transition point, 
“trans, Change very considerably. In series 4 a decrease in the value of the 
Reynolds number Re, at the transition point can be observed. 

It was pointed out before that the calculations of the pressure wave of 
series 2, using a coefficient of friction corresponding to Re, in steady 
developed flows, show that the assumed coefficient was too high in the 
initial part of the wave and too small at the end of it. Fig. 14 shows a 


comparison of the values of the coefficients of friction in non-steady flows 
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as calculated for series 2, with the steady flow values. The steady flow 
coefficients as obtained from Nikuradse’s formula are shown by straight 
lines as they are constant for a given Mach wave (characteristic). 
Fig. 14 confirms the conclusions drawn from Figs. 2 and 3 and indicates 
the fallacy of the application of steady flow coefficients to non-steady flow 
calculations. For laminar boundary layer the non-steady coefficient of 
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friction is shown to be smaller than the one corresponding to steady flow 
and based on Re;. One may expect, however, that it will not always be so 
because the influence of acceleration in non-steady flows, as indicated by 
Moore’s correlation equation, tends to increase the value of the coefficient 
of friction. 
4. Conclusions 

1. The initial boundary conditions have a very pronounced effect on the 
pressure wave produced. In the case investigated the cooling effect of the 
piston face and the tube walls and the effect of friction influenced the 
produced pressure wave to a large extent. 

2. The influence of friction in the range investigated was found to increase 
the pressure, temperature, and entropy of a non-steep compression wave 
produced by a moving piston above the values corresponding to an isen- 
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tropic wave, and to decrease the pressure, temperature, particle velocity 
und entropy of a moving wave. The decrease in entropy is due to the 
decrease in the lengths of the paths of particles with distance travelled by 
the wave 

3. In the case of air compressed by a piston through a pressure 
ratio of 1-6 in a tube which has a temperature equal to that of air before 
compression, the influence of heat transfer to the walls on the pressure of 
the wave produced is appreciable, but it is small in comparison with 
frictional effects when a moving wave is considered. 

t. The coefficient of friction in non-steady flows depends on the boundary 
layer formation and its value changes according to whether a laminar, 
transitional, or turbulent boundary layer has been attained between the 
cross-sections considered. It is, therefore, necessary to know the history 
of the motion in order to be able to predict the value of the coefficient of 
friction. 

5. The value of the coefficient of friction in a moving undisturbed wave 
lepends on the Reynolds number based on the distance travelled by the 
mean particle considered. This disproves the findings of Jenny (8) that the 
values of the coefficient corresponding to the Reynolds number based on 
diameter for steady developed flow conditions should be used. 

6. It appears that for a laminar boundary layer in a non-steep com- 
pression wave the coefficient of friction may be correlated using Moore’s 
correlation equation. This has been indicated by a comparison of correlated 
values with those following from Langhaar’s analysis of flow in the ‘inlet 
length’ toa tube. It follows that the influence of acceleration is to increase 
the value of the coefficient of friction. 

7. The values of the correlated coefficient of friction as presented on 
Fig. 13 are applicable to calculations of non-steep compression waves in 
compressible fluids moving in pipes under similar conditions (negligible 
heat transfer), in which a laminar or transitional boundary layer has been 
attained. In the case of a transitional boundary layer the coefficient may 
differ somewhat from the value shown on Fig. 13, depending on the value 
of the Reynolds number Re, at the transition point. Equations (21) and 
26) should be used to obtain the actual values of the local coefficients of 
friction. 

8. The calculations of a moving pressure wave for which equations have 
been derived in (1) and (2) need not be very laborious. Only a slight error 
ls introduced if a large net is used. 

9. The transition from the laminar to turbulent boundary layer was 
found to begin when a pressure ratio of about P/P, = 1-35 had been reached, 
P, denoting the initial atmospheric pressure. 
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10. As the actual flow is not one-dimensional, one must note the possi- 
bility of a radial flow due to the pressure loss in the boundary layer, and 
a consequent influence of the tube diameter on the pressure loss. 
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ON THE STABILITY OF SOLUTIONS OF CERTAIN 
DIFFERENTIAL EQUATIONS OF THE FOURTH 
ORDER 
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SUMMARY 
The Routh-Hurwitz criteria for the stability of solutions of linear differential 
uations of the fourth order are generalized for certain types of non-linear differen- 
tial equations of the fourth order by the use of Lyapunov’s function V. The method 
s similar to that of Barbasin and Simanov for third order equations, but yields 


somewhat less satisfactory results. 


1.1. Introduction 

Mvcu work has been done on non-linear differential equations of the second 
order, and it has been shown that the conditions for stability for linear 
equations of the second order can be generalized without much difficulty. 
For higher order equations the Routh—Hurwitz criteria for stability for 
linear equations with constant coefficients are well known, and also the 
results of Lyapunov which express criteria for stability in terms of a 
juadratic form V. Comparatively little has been done to extend these 
methods to non-linear equations of higher order, but Barbasin (1) and 
Simanov (2) have each applied a different type of Lyapunov function V 
to a third order equation obtaining different generalizations according to 
the function used. The object of this paper is to apply their methods to 
fourth order equations. Unfortunately the Lyapunov functions for fourth 
order equations are less easy to calculate and less manageable in generalizing, 
so that the results so far obtained are less satisfactory. 


1.2. History of the problem 
M. A. Aizerman (3) suggested the following problem: consider the system 


of differential equations 


lx 
= + Gut, 4 ae n—1 
dt =1 : 
, (1) 
AX 
: > a,,%;+f(2,), 1 k s n | 
k 
dt =a ss 
where a;; is constant, and suppose that if 
f(x,.) hzx;., r<h- B, (2) 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 2 (1956) 





186 M. L. CARTWRIGHT 


the roots of the characteristic equation have negative real parts so that 
the solution x, = x, = ... = 2, = 0 is stable. Does this remain true for 
more general f(x,.), provided that 
vty, <f(x;,) < Bar,? 

Erugin (4) and Malkin (5) dealt with the case n 2 by constructing a 
Lyapunov function for the system defined by (1) and (2) and then general- 
izing it. This method was used by Barbasin (1) and Simanov (2), but it leads 
to generalizations of a slightly different type from that suggested by Aizer- 
man, although Barbasin’s work practically covers one important case 
with nm = 3 and k = 1 or 2. The equations considered are all of types 
which occur in connexion with automatic control systems. 


1.3. Routh’s criterion 
All solutions of the differential equation 
(4) > ve —_ 
uA, €+A,%+0,%+0,2 = 0 (1) 
tend to the trivial solution z = 0 as t > 0, if and only if the roots of the 
characteristic equation have negative real parts. Routh (6) showed that 
this was the case if the first coefficient of the corresponding Sturm functions 
are each positive. That is to say if 
y 2 2 ; 9 
a,, Az =a,a,—a3, A, = a,a,a,—a2—aZa,, Aya, 4,4, (2) 
are all positive. This implies that 
9 
1, Ag, As, Ay, Ay (3) 


are all positive, and we assume henceforth that this is so. 


1.4. Lyapunov’s function 
The equation 1.3(1) may be written as a system 

z=%, @ 2 £=e we, ww —A, W—Ag2—Ag Y—Ag 2. (1) 
Lyapunov (7)+ showed that, if the real parts of the roots of the characteristic 
equation are all negative, corresponding to each positive definite quadratic 
form U(x,y,z,w) there corresponds a positive definite quadratic form 
V(x, y,2,w) such that O ae — i 
The method has been developed by Malkin (5) and others and holds good 
for positive semi-definite (7, and BarbaSin (1) and Simanov (2) use a 
function U which is a perfect square. This simplifies the calculations, and 
simplicity seems essential for effective generalizations if only to avoid 
being overwhelmed by formulae. 


+ Liapounoff, p. 277. 
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O that 2.1. BarbaSin’s V 
ue for . : 
‘or the equation . : . 
I | ita, #+(%)+f(x) 0 (1) 
Barbasin used 
| 2V 2+a, y)*?+20(y)+2f(x)y+ 2a, F(x), (2) 
ting a | . 
neral- | where F(a | f(~) dx, D(y) | d(y) dy, and # = y, ¥ z a before. 
| ] 0 0 
eads | : 
[his corresponds to the function 
\izel 
eine 2V, = (z+a,y)?+ 4, y?+ 2a, ry+a, a, 2* (3) 
types for the lineal equation 
t+a,%+A,%+AgxX 0 (4) 
obtained by taking U to be a multiple of y?. By skilful choice of functions 
> in (1) and (2), BarbaSin has obtained a most satisfactory generalization in 
which only the most natural restrictions are placed on ¢ and f, viz. V + oo 
ys at y?2+-22 = o, f’ continuous, and the three conditions 
t the $(y) 
} a 0, Nxeyie > Y, a J f'(x) > 0, 
that I 1 . 
, y 
tions ae 3 ° P . ° ~ 
} which obviously correspond to Routh’s criteria for (4) in the form 
; a, (), a, > 0, A, A, Ag—Azg > 0. 
2 In general we shall use V;, to denote a Lyapunov function for which 
U, where U is a multiple of y?, that is #*. This type of Lyapunov 
(3 function seems most suitable for considering generalizations of the coeffi- 


) cient a, in 1.3 (1). 


2.2. The fourth order linear J), 
Since a positive definite V;, is known to exist we may assume that it is 
given by 


|W, = by(w-+-Az t+ By+Cx)?-+ ky(e-+ Dy-+ Bx)?+ kg(y+ Fx) +hy2%, 
tic | Where ky, ky, ky, ky are positive, and that Vs is a constant multiplied by y?. 


Then by 1.4(1) 


rn 
l, k, (w+ Az+ By+ Cx){(a,—A)w-+ (a,— B)z+ (ag—C)y+a4 x} + 
k,(z+ Dy+ Ex)(w+ Dz+ Ey)+kg(y+ Fx)(z4+- Fy)+hy xy. 
; / Hence equating coefficients other than that of y* to zero, we have for the 
nd | coefficient of 
oid (i) w k,(a,—A 0, 
(ii) 22: k, A(a,— B)+k, D Q), 
k, Ca, 
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From (i) and (iii) A = a, and C = 0, and these facts will be used jn 
what follows. 


(iv) wz; —k,(a,—B)+k, = 0, ) 


(v) wy; —k,a,+k,D = 0, 

(vi) wx; —k,a,+k, FE = 0, 

(vii) zy; —k,{a,a3+B(a,—B)}+k(E+D2)+k, = 0, 
(vill) 2%; —k,a,a,+k, ED+k, F = 0, ) 


(ix) yx; —k, Ba, +k, B2-+k, F2+k, = 0. 


From (ii) and (iv) D = a,. From (v), (iv), and (vi) we have 





k k. a — ee 
: , B= a,——, BE — 14, 
ay, a3 ay, a3 


Then from (viii) F = 0, and by (vii) k, = k,A,/a?. Lastly, by (ix), 


ky = k,a,A,/(a,a,), so that putting k, = a}ag, we have 


\ 9 9 
a. p a,a = 
hy 4 aa,(: | ayy +t) ah 


2), = a2a,!w La, z-} (as. 
| ay 3 


+a, A,y?+a,a,A,2". (1) 


2.3. A generalization of J}, 

The effectiveness of Barbasin’s generalization is due to the fact that a, 
and a occur in 2.1 (3) in association with the corresponding variables y 
and x respectively so that he could replace 4a, x? by F(x), a,x by f(x), and 
5a, y” by M(y). In 2.2 (1) each coefficient occurs at least once in connexion 
with some variable other than that with which it occurs in 1.4(1), and 
each coefficient occurs in connexion with at least two different variables. 
This makes generalization difficult and the most effective one which I 
have been able to find consists in replacing a,x by f(x) in 1.4(1), and a, 
by f’(#) in 1.3 (2). 

We therefore consider the equation 


4 +@,%+a,¢+a,¢+f(x) = 0, (1) 
or rather the system 
t=Y, Y=2 t=u, & = —a,w—a,z—azy—f(z), 2) 
with 
; a3\ \? a ‘ 
2V;,(x, y, 2, w) aids) wta,z+ (22 ‘yy +afa,(z+a,y+—f(x)) + 
| a,)" } ag 
+a, A,(x)y?+ 2a,a, A, F(x)—ai f(x), (3) 
where A ,(x) = a, a,a,—az—a? f'(zx). (4) 
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THEOREM |. Suppose that (i) ay, ag, a, are all positive, (ii) f(0) = 0, and 


(x) > 0 and f"(x) is continuous for all values of x considered, (iii) 
F(x) = | pla) ax > OO 
0 
as |x| > ©, (iv) A(x) > 8 > 0 for all x considered. Then for every V, > 0 


7 : , ; . y r¢ 
there 18 a domain Dy of the x. Y, 2, W space defined by Ve(x, y, Z, w) — Vy.* 


If Poo. Yo: Zs Wy) 18 any point of Dy, then the solution of (2) through P, tends 


to the trivial aie a y =z= w= 0, of (2) as t +> @, provided that 


f"(x)y| < d/a, (5) 

107 all x, yen D,. 
It should be observed that, since f”(x) is continuous, the solution is 
uniquely defined by its initial conditions. Since f(0) = 0 and f’(x) > 0, 


ve have f(x)/a 0, and so 

. ix 3 f2/> 3 

2a, a, A, F(x)—ai f?(x) 2 | flx)faya, A 3— aif '(: 

0 
2 | a, f(x)Aq(x) dx > da, F(x) > 0 as |x| > 00. 
0 
Hence in (3) V, > 0 unless x = y= z=>w= 0, andV,>o as 
xt y?+-22+ w? + oo. 

It follows that for all ¥, > 0 the inequality V,;, < \, determines a domain. 


For any solution of (2) straightforward calculation gives 
Vy, = —ay as{A q(x) +40, yf "(x)}y? (6) 

D,, provided that (5) holds, and the rest follows more or less as in 
Barbasin’s paper. Let x(t), y(t), 2(t), w(t) be the solution through the 
point P, of the system (2). Then by (6) V,(t) Vfa(t), y(t), z(t), w(t)} is 
non-increasing and non-negative, and V,,(t) tends to a non-negative limit 
V,(00) as t>oo. If V,(00) 0, the solution of (2) satisfying the given 
initial conditions tends to the trivial solution z = y = z = w = Oast>o 
so that we have the required result. 

If V,,(00) > 0, then the surface V,(x,y,z,w) = V,(00) contains all the 
limiting points of {a(t), y(t), z(t), w(t)}. Let P, be a limiting point, then it is 
known that the sathitlinn through P, at ¢ = 0 lies in the surface 


V(x, y, z,w) = V,(00). 


Hence J 0 at all points of this solution, and this is only possible if 


either A,(x) a, f"(x)é or y= %= 0. Since f"(x) is continuous, it 
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has a maximum M on V, = V,(o), and so the first condition implies that 


for this solution 28 
r| > - ~ 
ai M 
If this is the case |x| > 00, which is impossible, and if y = 0, 2 is constant, 
and therefore x = 27, 40, y= y, = 0, 2 = 2, = 0, w = w, = 9, which 
contradicts (2). 
3.1. Simanov’s V 
For the equation 
+(x, £)He-4 a,i+a,x = 0 (1) 
Simanov used the function 
Vy = ia,z+a,y)?+4a,(a, 2+ a, y)?+3C(x, y), (2) 
eS | 
where J (x,y) 77 | u(x, y)y dy — a, ys 
0 
and assumed, in addition to the three natural conditions, a, > 0, a, > 0, 


W(x, y)ag—a, > 0, which correspond to Routh’s criterion, the condition 


y < 0. (3) 
Cx 


The function in (2) corresponds to the function 
2V5 = (ag2+-3y)?+a,(a,¢+4, y)?-+a, Asy? 
for the linear equation 2.1 (4), and here J; = —U, where U is a multiple 
of z*, although for the V, given in (2) there is, in addition, a term in V, of 
the form ¥ 
* Cus 
“Ay Az Y y dy 
Cx 


. 


which is non-negative when (3) holds. We shall use Vy, to denote a Lyapunov 
function obtained from that for a linear equation with U a multiple of 2’. 


3.2. The fourth order linear J, 
This time we suppose that 


) 


2Vz k,(w+ Az+ By+ Cx)?+k,(z+ Dy+ Ex)?+kg(z+ Fy)?+-kyy’, 


where k,, ky, kz, ky are positive, and try to choose the coefficients so that 
> 


k U is a multiple of z?. By 1.4(1) 
| 14 k,(w+ Az+ By+ Cx){(a,—A)w-+ (a,— B)z+ (ag—C)y+a,x}4 


+ko(z+ Dy+ Ex)(w+Dz+ Ey)+ks(z+ Fy)(w+ Fz)+khyyz, 
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and so we have A a,,C 0, as before, and 

(i) coefficient of y? k, Ba,+k, DE 0, 

(ii) coefficient of wz k,(a,.— B)+k,+k, = 0, 
(iii) coefficient of wy k,a,+-k, D+k, F 0 
(iv) coefficient of wa kya, ky E 0) 

(v) coefficient of zy k,{a,a,+ B(a,— B)}+k,( E+ D?)-4 

T ks i tk, 

(vi) coefficient of za ky, a, + k, DE 0, 
(vii) coefficient of yx k, Ba,+k, E? (0. 
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From (i) and (vi) we have that Ba, = a,a,, and then from (iv) and (vii) 


we have EF = B, and so a,k, = k,a3. From (i) and (iv) D = a,, and so 
from (iii) it follows that F 0. Hence we have 
| kay Q,4, 3} hy Ay ky pas 
| As a,) A, 4 az 
and, putting /, a,az, we have 
a,a : a,a : 
SD) i 2 a 4 i, 3] + _] a "6 oe »2_} 2 
2, = a, asl w+a,z : /) a3(z %4y+——2) 4 a, A,z*+a,a,A,y’. 
3 3 / 


3.3. A generalization of V. 


\cain each coefficient occurs a most inconvenient number of times, and 


is associated with variables with which it has no obvious connexion. 


However, a, occurs only twice, each time in A, in connexion with y or 


so we replace ta,y*? by | yu&(y) dy = V(y), and in other cases replace a, 


t] 
by u(y). We therefore consider the equation 
{ 


ca uh( Y)a aaX--a 


. ) 
3 fl 0, 


or rather the system 


y j u u a, w—w(y)z—dg y—a,a 
with 
a,a - a,a 2 
QV. ry 1,43) U 1 l ty a3(2 ayy . ty 
As As 


a.(a,deu(y)—az—a 


ee 3 l 





2 42 D2 Wy 2 2 2 
a,)z*+ 2aja,a,V¥(y)—a,a,(a3+ajza,)y*. 


2) 


(3) 
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THEOREM 2. Suppose that (i) a,, 43, a, are all positive, (ii) y(y) 18 positive 
and us'(y) continuous, and (iii) 


A,(y) = a,a,%(y)—a3—aja, >8 > 0 (4) 


for all values of y considered. Then for every V, > 0 there is a domain D, of 


the x, y, z, w space defined by pi z,w) <V. Lf Po(%. Yo: 2 Wo) 18 any 
point of Dp. then the 90 of ( sonia P, tends to the trivial solution 
x=y=z2=w= 0, of (2) as t+, provided that 

as\p'(y)| |z| <8 (5) 


for all y, z in D,. 
Since (4) holds 

y 

, 9 ¢ 
2a} a,a, V(y)—a, a,(a3+aja,)y” = 2a,a, |A (yy dy 
0 
; 7 
> A, 4, dy? > ©, 


9; 9 9 


as yoo, and so Vy is non-negative and V > 0 as x7+-y?+-27+-w* > ~, 
Hence, as before, for all VY, > 0 the inequality Vz < Vj, determines a 
domain D,. 
For any solution of (2) straightforward calculation gives 
Vs = —a,a,{A,(y)— sat y)ziz27 < 
so long as (5) holds. As in Theorem 1 V.{2(t), y(t), z(t), w(t)} tends to a non- 
negative limit V,(00), and if V(co) > Oa oa must lie in the surface 
V, = V,(00). This can only be the case if either z = 0, or 
A,(y)—ha,y'(y)z = 0. (6) 
If z = 0, then y is constant, and either 2 > +o, or y = 0, and we have 
a contradiction as before. On the other hand, if (6) holds it implies that 
2z| = |y| > 28/(a, 0), 
where M is the maximum of #’(y) in D). If this is the case |y| > 0 which 
is impossible for V, = — < Vy. Hence V,(00) = 0, and the solution of 
(2) through (2p, Yo, 2, @») tends to the origin which is equivalent to the 
result stated. 


4. Comparison of ;, and J, 

The different effect of using V, instead of V;, when both are possible may 
be seen quite clearly in the case of third order equations by attempting to 
use V;, for Simanov’s generalization. For simplicity we may suppose that 
#% is a function of # only. Then keeping 
by #(y) in 2.1(3) we have 


ay, 4, constant and replacing 4, 


B= —(d(y)ag—ay)y? +a, b' (y)a?z + {z+b(y)y eh’ (y)yz 
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If, corresponding to Routh’s criterion, we have My, d, positive and 


Yjas e > oO > QO, 


then JV), 0 when 
b(y)|\dsx°z+-tz+ply)yjzy| < dy? 
for the values considered, and this is not necessarily satisfied unless 
wb (y) = Oly?) 


as y > 0, but no such condition is required in Simanov’s proof. 

On the other hand, by a skilful choice of Vg it is possible to handle a 
function % of x and y only, when the dominant part of Vz is a multiple of 
‘although since ¢s(x, y) is not a function of z no condition can be imposed 


on to make it tend to 0 as the dominant part of V, tends to 0. 


5. Other seneralizations 
The method seems suitable for application to a fourth order equation 
which any coefficient is a function of 2, #, #, or % as the case may be. 

For by taking U’ for the linear equation to be a multiple of the corresponding 

variable 2°, y®, 2°, or w®, whatever the case may be, we can obtain a 

Lyapunov function V such that the dominant term of V is not greater 

than —U’, and we can impose a condition which will make V < 0. For 

example, if a, is replaced by a function y(x), we should have to construct 

’ such that the corresponding linear U is of the form 8z?, and then the 

general V would involve a term similar to —U together with terms in 
r)y. Hence y'(a) would probably have to satisfy a condition of the form 

O(x?) as x + 0. In many systems the function involved is almost 
exactly linear near the origin, and so this would not be a serious drawback, 
but the fact that for third order equations a different V gave a result 
independent of such a condition makes one hesitate to pursue this line. 

The calculation of Ve and JV. 


, for linear equations seems to be somewhat 
simpler than those for more general V involving two or more variables, 


though this might be more advantageous in some respects. 
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90 SUMMARY 
— The problem studied is that of the vibrations of a flexible beam supporting a 
7 moving locomotive; associated with the locomotive is an alternating load caused 
the ‘hammer-blow’ effect. A fundamental differential equation is established 





und an approximate solution of this equation is derived by making use of the 
Fourier expansions of certain Jacobian elliptie functions. The predicted ‘dynamic’ 


{ deflexions of a bridge of 80 metres span are in good agreement with measured values. 


critical speed of the locomotive, for which the dynamic deflexion reaches a 
| maximum value, is predicted accurately. 


1, Introduction 

VIBRATIONS in railway bridges are commonly caused by the ‘hammer- 

? blows’ associated with steam locomotives, in which counterweights are 
attached to the driving wheels to minimize the horizontal inertia effects 
of the moving piston. Unfortunately the counterweights may give rise 
ilso to a vertical alternating inertia force, which is transmitted to the 
deck of the bridge. When the locomotive crosses the bridge at constant 
speed, the hammer-blow varies in a simple-harmonic fashion, the period 
being the time required for the driving wheels to complete one revolution. 
In discussing the behaviour of the bridge under the hammer-blow we 
assume that the locomotive is a concentrated mass, travelling at constant 
speed along a uniform beam which is simply-supported at each end. Thus 


no account is taken of the ‘springing’ of the locomotive on its wheels. 


2. Fundamental differential equation 

A uniform horizontal beam of length /, flexural stiffness EJ, and total 
mass M, is simply-supported at each end (Fig. 1). Consider the vibrational 
effects caused by a locomotive of mass m, travelling along the beam at 
constant speed v. The vertical force due to the hammer-blow at time ¢ is 

H sin(wt+e), 

where H is the maximum force, w is the angular velocity of the driving 
wheels, and H sine is the force at time t = 0. 


This paper is a summarized version of ‘The forced vibration of bridges under the 
g loads’, by the late K. Mise and 8S. Kunii (Kyoto University). The 
summary has been prepared by A. H. Chilver (University of Cambridge). 


pF action of movin 





[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 2 (1956)] 
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Suppose that y denotes the downward vertical deflexion of the beam 
at a distance x from one support at time ¢ and that the coordinates of 
the locomotive at any time are (X,Y). We may take the deflexion of the 
beam to be of the form 











dl 


Y 
Frc. 1. 


in which the /,(¢)’s are functions of ¢ only. By evaluating the total kinetic 
and potential energies of the system, we find that the Lagrangian equations 





of motion give the differential equations 


d*f. ; J? , 
=e +. 2h = + (ra)*f, = y|g- Hsin (ot +e) a >, f sin sBt |sin rft, 





(1) | 
where r takes the values 1, 2...., 20, h is the coefficient of velocity damping, 
and «a, 8, y are given by 

4h 
R mEI 


¥ MB 5 p (2) 


~ 
~~ 
. 
| 
\ 
— 
i 


3. Approximate form of the differential equation 

In comparison with f, it can be shown that the value of /. is of order r~‘ 
so that the values of f,, fs,... are negligible: moreover yf? sin*St/a* is also 
negligible compared with unity. An approximate form of equation (1) is 
then 


(1+ysin?8t) a’; 2(h- yB sin Bt cos Bt) qf, xf, 


dt? | re | 
wl 


Ao+; saad c) sin Bi (3 
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Let f; n, Bt = €, and write 
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eam | Then equation (3) becomes 
1S yt | d*y ro We dy a. we ¢ 7 
Fthe | — 4 [2A0Q(é)+ P(é)| et AVO(E)n = F(ENQ(E), (5) 
| Ge as 
| in which 2y sin € cos €& 
P(§ ) =o ss aay (6) 
l+ysin°é 
l s 
| V(E&) me 9 (7) 
l+ysin*g 
F(é) = |ygo +2) sin(ué +) |sin €. (8) 


The function » gives the deflexion of the beam at mid-span at any instant f. 
Since the bridge is initially at rest the initial conditions are given by 
dy 
dé 


0, att 0. 


[he equation for free vibrations is obtained by equating to zero the right- 
ind side of equation (3); the lowest natural frequency of the beam is 
neti riven by 9 ° 
4 ; h?\3 
tions Ne l )*. 
2a “i 
In practice the value of h?/a* is small in comparison with unity, and we 


; ; ; a ie ~ 
, may write with sufficient accuracy a* 277No- 


|) | 4, Solution of the approximate equation 


ping, | Consider firstly the solution of the equation 


dy di 
rag + (2A0Q+-P) et On 0, (9) 


btained by equating to zero the right-hand side of equation (5). On 


making the substitution 
wp gad 7 wexp| | (AoQ+ 1P) dé], 


equation (9) becomes 





124 ul 20 \oQ +1 PY? (A dQ l at 0. 
a, dé 2dé 
If\ is large, an approximate solution is 
u (Q—oa?Q?) ‘exp|ia ( (Q—oQ?)! dé, 
1), p. 490). But since |Q l,o*? <1,and PQ dQ/dé, a sufficiently 
ccurate value of 7 1s 


Q' exp| —Ao { Q délexp|ia | q dé. 
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Two independent values of 7 satisfying equation (9) are 


n = @! exp|—Ao | Q dé|cos|) | Q! dé| 


ne Qi exp| do { Q dé|sin|A | Q! dé| | 
The solution of equation (5) satisfying the initial conditions is 


é g 
mC) F(C)Q(C) dC »s(é) 
. ni(C)no(Z)—4(S)o(Z) 


)Q(L) dv 
n = n2(€) PUC)O(G) dt 


Nol 
mit 


5. Undamped forced vibrations 


equation (5) reduces to 


d*y dy , " 
oe 4 t A¢ k . 
de dé v7 Q 
é 
Write @M(€) = (1+ ysin7€)-+, V(€) = | (l+-ysin®Z)+ dé. 
0 
Then equations (10) become 
m = O(é)cos[AF(E)], 


and equation (11) may be rewritten 


O(€)sin{ AV (E)], 


tw 


I 
” 1 PENI I,+2H,1,), 


where  & [ @(C)sin ¢ sin| A‘’(€)—A'F'(Z) | dZ, 
0 


‘4 
$ 





and ly ] (C)sin ¢ sin| AY (€)—AV(C) |sin(ul+-e) dé. 
0 
_}"3 

If we write k’2 ; * q a am 

] 

y 
A(1—2q)? 

and 6 = 7X /l, 


accuracy, by the relations 


i = ,(sin 4 -V sin ; : 


)n aa ~n(L)9(f) 


(10) 


(11) 


When there is no damping present, h = o = 0, and the differential 


(13) 


(14) 


(15) 


(16) 


(17) 


(18) 


then it is shown in the Appendix that J, and J, are given, with sufficient 


(19) 
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| y 
I, S 1 (da, bi G(u 2) 1) -G(uw—2r +] )} t 
2 aa 
(a,—b,){G(u+2r—1)—G@(u+2r+1)}], (20) 
I which 
— cama . cee - 
G(N) (1—N2V2)-1| cos(N@+e)+ NV sin esin 7 — 008 € COs 


where V takes the values (u+2r-+-1), and 
































a, l uq?a, (1 uq?)a, | 1 uq?a}. ‘s [7 om Ot 5], (22) 
bo U b. Lug?b) »+(1—3yq")b.+}u9qb,.. [r= 1 to 5}. 
(23) 
In addition, we have 
Lo 2 0O 12 0 ] (PH) g? ay — 
nm i 4 Q 15 (74 \q3 a; 
} 0 0 3 O (*H)q4 a - 
0 Lg 0.2 (RH) gS a. ma} 
2 1 32 4/9 3 
00} 0-30 (“E)q8 a4 
1 2 9 , 
0 8 @ | 0 33 ] (eq : | as 
1d 0 3 0 2 0 35 ("Hg ] bi —2yq 
1 3 2 2 ‘ ’ 
50 ¢ YU 5 0 ( *)q° bs 
l 21 2 , oR 
0 + 0 0 2 ( H\g4 b; (25) 
0 O 0 & 0 (7#)q® bi 
; 
0 O O 0 a (*#)g® b; 7 
Lg" 








6. Damped forced vibrations 
When damping is present the values of 7, and yn, are given with sufficient 


accuracy by the relations 


Y 


1 O(é)exp[ —A, AV (E) |eos[ AV (E) | | 


(26) 
no = O(é)exp[ —h, AF (E)]sin[AV(E)] J 


where hy ma/2K, (27) 


in which K is the complete elliptic integral of the first kind of modulus k, 


where 
F 


= (28) 


k2 = 
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Equation (11) gives 


l " ~ , , 
heads O(E)exp[ —hy AV (E) (yoo 14 +2, 15), 29) 
, ] pall . 
where i= 5 LexP (lta 8 V)sin @—V sin(@/V)], (30 
and 
» eS... 
f= gy D, [er NGu—2r—1)— 62 +1)} 


+ (a,—b, KG" (u+-2r—1)—G" (w+ 2r+1)}]; (31) 
a,, 6, ave defined as before by equations (22) and (23), and 
G'(N) = {1—N?V?)[exp(hy 4/V )cos(.N6+-¢)—cos € cos(@/V)|4 
_NV(1—N?V?)sin e sin(@/V)+- 
+ 2NVhlexp(hy 4) V)sin(N6+-)—sin e cos(4/V)] 
—h,(1+.N?V?)cos € sin(6/V)} > 
{(1—N2V2)?4+- 2h2(1+ N2V2)}-1, (32) 


where takes the values (u+2r+-1). 


7. Comparison with experimental results 

The vibrations of a railway bridge have been studied by Sir Charles Inglis 
(2), who measured the deflexions at the mid-span of the Newark Dyke 
Bridge for a single locomotive travelling over the bridge at various speeds. 
The measured vertical deflexions at mid-span for four speeds of the loco- 
motive are shown in Fig. 2. Except for the damping coefficient, all the 
data required for a theoretical study of the vibrations are given by Inglis; 
the coefficient of velocity damping is found by comparing two successive 
amplitudes of free vibration occurring after the passage of the locomotive 
over the bridge. 

The relevant constants for the Newark Dyke Bridge and the 0-8-0 type 
locomotive are 


/ = 80 metres, M 467 tons, N = 2-88 cycles/sec., 
h/n, = 0-0925, «= 0, m 109 tons, 
H 0-315v? metre tons /sec., 
circumference of driving wheels 4-27 metres. 


Then y = 2m/M = 0-467, q = 0-0239, b= 37-5. 
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The values of (a,+-6,) and (a,—b,) are 


dg +by) = 0-342 (ay—bo) = 0-342, 
a, b,) 1-17 (a, b,) 1-15, 
a,+b,) = 0-604, (a,—b,) = 0-618, 
a,+be) 0-185. (a,—bs) 0-210, 
(a, b,) 0-0393. (a, by) 00-0519, 
(a b.) 0-00658. (a- b;) 0-O107. 
ans 4m.sec v= 24sec 
sec v= 2-40sec 7 
» DAANDA J 
= = ~ +} \ ae, 4 
ec v = 2-69sec 
- nd LAF 
m.se v = 6-20sec 
4 5 ¢ 70 a 9C 00 m 
Fic. 2 
[he values corresponding to 7 t andy = 5 are small, and may be safely 
neglected. We have further that 
v 417) metres/sec., ho 0-0133. 
For the lowest speed of the locomotive (v 9-14 metres sec.) the deflexion 


t the centre of the bridge 1S 


mm 14-0 sin 6— 0-5 cos 38:58-+-0-7 cos 40-58 — 0-5 cos 42°56- 


0-3 cos 44°:50-+ 0-2 sin 44°50 0-4 exp| 0-6076|sin 45-68. 


The first term on the right-hand side represents the ‘static’ deflexion of 


the bridge 


for other speeds of the locomotive the deflexion equations take 
$1! i] ir forms 
The theoretical deflexions at mid-span for four speeds of the locomotive 


re plotted in Fig. 3. The theoretical and experimental curves (Figs. 


ind 3) compare favourably. At the lowest speed of the locomotive the 


maximum amplitude of ‘dynamic’ deflexion of the bridge is small, and 
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occurs just after the locomotive passes the centre of the bridge. At the 
next speed the maximum amplitude is large, and is nearly half the static 
deflexion; the maximum value occurs when the locomotive passes a point 
two-thirds of the distance along the span. At the third speed the maximum 
amplitude of ‘dynamic’ deflexion is still large, and occurs just before the 
locomotive passes the centre of the bridge. At the highest speed the 
measured ‘dynamic’ deflexion varies erratically, but the maximum ampli- 
tude is small. 


eae v= 9-l4m.sec v= 2-14sec 
v=10-2m.sec v= 2.4Osec — 
0 es A wit 4 
v = 1|.5m.sec™' v = 2.69sec"' ; 
0 in APA 
fea cea — = 
| iT ET AE OE | 
200 , ; , fcanes 4 1 4 — + 4 J 
_ v=26.5m.sec v = 6.20sec"' a 
—— donee eee eeccuinctesio’ 
100 | ee ees _— = ian sae ——— 
200° 4 4 —_4—___— + 4 4 — nn) 
0 10 20 30 40 50 60 70 80 90 10Om 
Fic. 3 


8. Dynamic deflexion 


The total deflexion at mid-span is composed of the ‘static’ deflexion ng, 
and the ‘dynamic’ deflexion 7, and we write 


"7 IsT 1p: 
From the theoretical curves of Fig. 3 we see that 7, attains its maximum 
value when the speed of the locomotive is about 10-2 metres/sec. To find 
this critical speed more precisely it is necessary to compute the dynamic 
deflexions for speeds between 10-2 metres/sec., and 11-5 metres/sec. The 
results of a number of calculations are as follows: 


Loco. speed v (m./sec.) 9:14 10:2 10:3) 10:5 106 108 11-5 26-5 
Maximum 7p (mm.) 2-0 6-7 V2 8-6 9-0 Yo 4:7 1-1 
(Max. 7p)/(Max. 73) 14% 48% 51% 61% 64% 55% 34% 89 


From this range of values of », we see that the critical speed is about 
10-6 metres/sec., and that the maximum dynamic deflexion at this speed 
is about 64 per cent. of the maximum static deflexion. 
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9, Comparisons with Inglis’s analysis 


Sir Charles Inglis (2) considered the deflected form of the bridge to be 


7X 
? ‘(t)sir 
y = f(t)sin j 


} where f(t) is a function of ¢ only. For damped forced vibrations he estab- 


lished the differential equation 


[4 Ff ooTanet a ame. Quvtl df 
I El ? ( eos = po[ hme sin ut if 
| M MP? l }y\ Ml 1 jdt 


| 2not\| 2f 2H. aw 
I m\'—°° r)| ae yp sinetsin 7) _ 


For the particular case of no damping (kh = 0) Inglis took the forced 
oscillation to be 
— rv 
f(t) S A, c0s(w-4 It 
This relation does not satisfy the initial condition that f = 0 when t = 0. 
\ free oscillation of the form 
— rm 
I(t) = B,005( 1 ary \e 


is superposed to correct for this discrepancy; w, and B, are determined by 
trial-and-error methods. Inglis takes terms up to r = 11 in the first of the 
above series, and up to r = 6 in the second. 


APPENDIX 
Evaluation of J, and /, 


We first introduce the quantities 
k2 ) : k’2 l k2 7 (34) 
Then we may writ¢ 


W(é { (1+ysin2Z)-# dz k’\ K 4 f a k* sin’) do}, 
where A is the complete elliptic integral of the first kind of modulus k. In terms of 


the Jacobian am-function, we have 


l 
lar am| 1 V6) K|. (35) 
Now Il, @(f)sin Z sin[A'’(€)—AV(f)] dé. 
Put us pt (¢), (36) 


where p w/2Kk’. (37) 
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The upper and lower limits of the integral J, are 
y= pT (0): = 0; xs = p¥(é) = @ (say). 


Then the transformation defined by equation (36) becomes 


2K 
C= thr4 am| ~(p- 17). (38) 
Suppose V pla; (39) 


then in terms of Jacobian dn and en functions, J, may be written 


 F me | an PE(y sn)en| =" (y 1x) |sin|F (0 | a. 


0 


But the en and the dn functions may be expanded as Fourier series of the forms 


en[ > (y 7) at qsin 35 ™ 
TT n Kk l1+q ]1+q 


9K | - Dien ne § 2a I, 
an| = (4 in) 2 i goes , q- cos “4 m 
1+q* l+q 


where q exp[—7K’'/K], and is the constant frequently used in the theory of 


elliptic functions; g may be evaluated with sufficient accuracy from the relation 


q > 


((3), p. 486). In our problem the value of y, which is twice the ratio of the masses 
of the locomotive and the bridge, is known; k’ is then found from equation (34), so 
that g is a known constant. When y = 0,q = 0, and when y 0, q }; so that all 
relevant values of q lie in the range 0 < q }, but in practice y takes values which 
make g small compared with unity. Then J, reduces to 


rs) 


1 (2q\i(7\3 - Fo , 
B ok =) | oy | (sin ys — 2q sin 35 ..)sinl E ( us ] dis. 
The quantities kh, k’, K, q are related by the equations 
(2kK /77) 2q*(1+q?+q°+...), 
(2h’ Kx /7r)? (1—2q+2gq!-+...) 
1/2q\3(7\3 
(3), ». £79). So that 7 2 l, 
stl iy) (x) 
1, l 
and I, | sinypsin| -(@ | a, 
p. | 
CS F ; . 
which reduces to iD V[sin @— Vsin(@)) (1—V?). 
p 


The value of p is given with sufficient accuracy by 


so that V p A 


; 
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} 
The function V’ is linearly proportional to the speed v of the locomotive, and is given 
pproximately by 1/A. The final approximate form of J, is 
l , ‘ 
I [sin @— V sin(6/V’)]. 
A 
38 It remains to discuss the function @ in this relation; the am function may be expressed 


: Te , | q’ sina Z/K 
Z 1) tan i q?' leos7Z K 


2h howd 
4), p. 285 So that equation (38) may be written 
c w—2tan 1(q sin 2y }q*sin 45), (40) 
’ tothe required degree of accuracy. Now 
6 pt (&) (us)¢_¢ é, 
approximately, J = ¢ Bt = wX/l. 
fact, 8 indicates the horizontal position of the locomotive on the bridge. Now 


5 I } D(C sin € sin[A4(&) AN(f) |sin(ul e) dl. 


As before, apply the transformation given by (38). Then 


l “ 
| sinpsin| + (6 J) sinquz e) dus. (41) 


n(ul +e in(pys + €)cos[ 2 tan—"(q sin 2x — $q? sin 4x5) | 
cos(pys + €)sin[2u tan—!(q sin 2ys— $q? sin 4y5)). 


s may be expressed with sufficient accuracy in terms of multiple angles in the 


> [(a b,)sin{(— 2r)b+-e}+(a,- b,)sin{(u- 2r)ys+-e}], 


where a,, 6, are defined by equations (22) and (23). On substituting this value of 


u¢+-e) into equation (41) and integrating, we get finally equation (20). 


Evaluation of J; and / 


When damping is present the functions 7, and 7, involve an additional exponential 


term. If 
u 
C lot amr ——K }, 
t= inten) 
é ¥(¢ 
u : 
then (1+ysin?Z)-! dZ k’ | nd{ >—K) du. 
0 0 


Introducing the Fourier expansion of the nd function, this integral reduces to 


7 “=qk 7 


¥(é) s Fe W@]+-. -¥(é), 
aK tia" lke *@) aK * (8) 
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since k’ < 1. The values of 7, and », are then given by equations (26); the integrals 
I, and Ij reduce to 


\ 


I, = [ ®(L)exp[hyAV(Q]sin £sin[AV(E)—AF(O)] aZ, 





0 
, ae eeiey, . By A. 
I, = | O(Ljexp[hy AV (Z)]sin f sinfAV(é) AY’(f)Jsin(ul +e) df. 
0 
As before, put p¥(¢). Then 
exp[—hpAT(f)] = exp[—hyts/V], 
; This 
z ] i ae - ] } for the 
and I; : | exp[hyys/I sing sin| + (6 | dis, may be 
0 ) of Sir C 
which is given with sufficient accuracy by equation (30). After transformation 
’ benann | 1. In} 
| becomes ff] | 
,» af , Z {. ., oseillat 
i exp[hy/V Jsin ys sin p (9 ys) |sin(uf +e) db, | 
} t 
Ps | 
which is given with sufficient accuracy by equation (31). 
; ) 
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BRIDGE VIBRATIONS 

By A. H. CHILVER (Engineering Department, Cambridge University) 

Received 21 March 1955] 


SUMMARY 


[his note describes briefly a number of points in relation to the Mise—Kunii theory 











r the oscillations of a railway bridge under the action of a moving load; these 
be of some interest in comparing this new approach with the earlier theories 
» of Sir Charles Inglis. 
atior ° 
1. ly his study of bridge vibrations Sir Charles Inglis established that the 
illations of a railway bridge are dominated almost entirely by the 
° Ty 
Psin 2aN é 
} 
Dv 
? —-——s s 
4 4B 
+ -— - vt ——? 
v 
Y 
Fig. 1. 


‘hammer-blow’ effects of steam locomotives (1); he found that for loads 
of constant magnitude, moving at typical speeds, the ‘dynamic’ deflexions 
ofarailway bridge due to oscillation are not large; for a uniformly distributed 
advancing load the deflexion is almost free of oscillation, and may be taken 
as the ‘static crawl deflexion’. These conclusions indicate that a more 
critical condition may arise when a single concentrated mass traverses a 
bridge; a theoretical analysis, to be of practical value, should take account 
of the ‘hammer-blow’ effects of the locomotive and damping effects in 
the bridge. This is the problem studied by Mise and Kunii (2). 

2. The paper by Mise and Kunii is essentially a mathematically more 
precise treatment of the problem studied by Inglis, who considers a 
simply-supported uniform beam of flexural stiffness EJ, total mass M,,, 
ind span /, Fig. 1 (the notation used is that of Inglis). The locomotive, 
of mass WV, crosses the bridge at constant speed v, and at an instant t¢ is 
a distance vf from one support. The force of the ‘hammer-biow’ at time tf 
has the value _° sin 27.NVt, where 27N is the angular velocity of the driving 
wheels. If «a is the downward acceleration of M at time ¢t, the total 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 2 (1956)] 
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concentrated load on the beam at that instant is (Psin27Nt— My). | damping 
omitting gravity forces. Inglis considers this concentrated force to he | value of. 





replaced by an equivalent distributed load of linear intensity cients of 
2 : a Sg OE, . 27x. . 38k . one of fo 
j (P sin 27Nt— Mx«){ sin j sin 27nt-+-sin - sin 47nt-+-sin j sin 6znt- wh Mise a 

nize the 
where n v/2/. Mise and Kunii point out that this series is not convergent: 


Inglis considers the concentrated load to be defined with sufficient accuracy | 
by the first term of the Fourier series, viz. 
») 7 
~ (P sin 27Nt— M«)sin 2rnt sin —, 
l l » consists 
and then argues that the primary and dominating component of the state | functior 


of oscillation is given by the equation 


(the not 


evaluat: 

oO" ran) raked) 2 : ; ; re 

EI°4 4 4nn, m oy m J (P sin 27Nt— M«x)sin 27nt sin 
ext ot ot? l l 


; 


where y is the downward vertical deflexion at any point at time ¢, , is the j 
coefficient of velocity damping, and m is the mass of the beam per unit 
length. On putting 472n5m = 7*EI/I*, where n, is the fundamental natural 
frequency of oscillation of the unloaded girder, and m/ = M,,, the total 
mass of the girder, and assuming a solution of the form y = f(f)sinza|], | in whit 


we have time f, 
ol. 22M 1. aM df 
ee tent |f : te] m4 — 5 
M d*f | ; 
+ |1+—— (1—cos 4rnt) |— | y being 
| u,! cos 477n No hme 
P . ‘ may s 
ly [cos 27(.NV —n)t—cos 27(N +-n)t]. unity, 
ies (84, 


The comparable equation of Mise and Kunii is derived from the Lagrangian } 7, 

J 
equations of motion, and gives essentially the same result. Mise and Kuni 
assume that the first term is given with sufficient accuracy by 47?njf. So 


although Inglis’s analysis employs a divergent Fourier series, the final | 
result is in agreement with that of Mise and Kunii. which 
3. It is in the solution of the differential equation that the two analyses , In pre 


diverge; Inglis takes the particular integral of the differential equation in | 
the form oer 
f(t) > [A,cos 27q,t+ B,sin 279, t], 

r x 


where g, = N-+-rn, and r is any odd integer, positive or negative. Inglis This i 
considers this to correspond to the forced oscillations of the beam; the free } More 


‘4 


° e ° . . - Ww ie 
oscillations are not important because they are suppressed rapidly by hict 
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Vy) |damping forces. The values of A, and B, are found by substituting this 


o be | value of f(t) into the differential equation, and then comparing the coeffi- 





ients of cos 27q,¢ and sin 27q,t on each side of the equation. In discussing 
ne of four numerical examples Inglis takes values of r from —7 to +13. 
| Mise and Kunii consider an interesting alternative approach; they recog- 
/ nize the use of elliptic functions in dealing with the differential equation 
rac d2y7 dy 


7 (+ [2A0Q(é)+ P(E)] 


pet OEn = FOE) 


the notation is that of Mise and Kunii (2)). Essentially their method 
onsists in deriving the solution of the equation from the complementary 


tate | function. In the case of no damping they reduce the problem to the 


valuation of the integrals 


( (C)sin f sin|AY'(€)—AVF (CZ) | dZ, 


l 
, 0 
the 
the } 
ut . 
” I, | O()sin ¢ sin| A‘’(€)—AV(C) |sin(ul+e) dé, 
ota 0 
“xl, | in which € defines the position of the locomotive on the bridge at any 
time t, and where 
g 
M(é) (1+-ysin*&)-¢, Y'(€) (1+-ysin2Z)- dZ, 
0 
y being twice the ratio of the masses of the locomotive and the beam. We 
nay solve J, approximately in the following way: y is usually less than 
nity, and when equal to unity varies only from 1-* to 2-+, i.e. from 1 to 
0:84. Let us suppose that ®(€) = 1. Moreover, suppose that Y(é) = &. 
ral : 
° Then J, reduces to 
ul a a n * - 
g I, | sin Zsin(A€—AZL) dZ, 
j U0 
vhich gives I, \ sin €—sin A€)/(A2— 1). 
ses, In practice A > 1; so that an approximate form of J, is 
nl 
\ a a » 
I, sin €— —~ sin A&}. 
é 
lis | This is essentially the same as the result derived by Mise and Kunii; their 
free | Morerigorous treatment of this integral involves the use of elliptic functions, 
by | Which is justified because J, can be evaluated similarly. Unfortunately it 
y | 2 


P 
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is impossible to evaluate J, accurately by a method as simple as that 
outlined above for J,; we have 

g 

I, = } | ®(Q)sin {[eosfr¥(€)—e—A¥ (0) — no} - 
0 
—cos{rv¥'(€)+-«—AF (CZ) +E} ] dé. 

An inherent difficulty in the evaluation of this integral is that AY’(2) and 
uf are of the same order of magnitude, and any approximations to Y({) 


| 
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Fig. 2. 


must take account of this; the way in which Mise and Kunii overcome this 
difficulty is of considerable interest. 

4. It is possible that some of the analysis in Mise and Kunii’s paper may 
be presented more simply. The important feature, however, is that their 
final result is probably simpler in its application to any particular problem 
than the method put forward by Inglis. The Mise~Kunii theory gives 
results which are in good agreement with the measured vibrations of the 
Newark Dyke Bridge over a wide range of locomotive speeds; Inglis found 
good agreement between his estimated and measured deflexions of the 
bridge for the case in which the locomotive was travelling at the ‘critical’ 


speed; at more remote speeds there are small discrepancies, which may be 
unimportant. The maximum ‘dynamic’ components of the central deflexion 
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f the Newark Dyke Bridge, caused by the passage of a single locomotive, 
we given in Fig. 2 for a range of locomotive speeds. The Mise—Kunii theory 
rives A Maximum dynamic deflexion about 30 per cent. greater than the 
ue estimated by Inglis. At the ‘critical’ speed, corresponding to the 
greatest value of the maximum deflexion, the deflexion is estimated by 
Mise and Kunii to be 64 per cent. of the maximum static deflexion. An 
ipproximation suggested by Inglis is that the maximum ‘dynamic’ 
leflexions of the bridge occur when the frequency of the hammer-blows 
sequal to the natural frequency of the bridge with the mass of the loco- 
motive concentrated at mid-span. This gives for the critical speed of the 
comotive ) Ny S| M, (M,, | 2M)}}, 

where %) is the unloaded natural frequency of the bridge, and s is the 
reumference of the driving wheels. For the Newark Dyke Bridge this 
pproximation gives v = 33-5 ft./sec.; the value predicted by Mise and 


Kunii is 34:8 ft./see 
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AN ANALYSIS OF TRANSISTOR PERFORMANCE AS A 
YTIATS UY T , ph . ) . > I wh ‘ 
FUNCTION OF FREQUENCY AND FOR REALISTIC 
GEOMETRIES 
By J. H. OWEN HARRIES (Warwick— Bermuda) 
[Received 9 June 1955] 
SUMMARY 

The boundary shapes of actual fused impurity transistors are not such that 
analytical solutions of the governing differential equations can be found for the 
flow of the carriers. A relaxation method is described for solving these differential 
equations for realistic boundary shapes and as a function of frequency. The linear 
small signal theory of transistor operation is used, 
1. Introduction 
EXIsTING mathematical analyses of transistors are in general limited to 
consideration of simple idealized shapes, such as the space between infinite 
parallel planes or parallel sided rods. The geometry of actual fused impurity 
transistors (Fig. 1) is entirely different from these simple idealized shapes 


EMITTER 












os BASE 
COLLECTOR 


Fic. 1. The geometry of a p—n—p fused impurity transistor. The base is 
n-type germanium. The emitter and collector are made of indium. 





and major factors in the design of these and other transistors depend 
critically upon details of geometry which cannot be taken into account by 
idealized shapes. This paper describes a method of applying relaxation to 
the study of transistors so that they can be analysed for realistic boundaries 
and geometries. 
2. Linear low-level theory 

The well-known linear small signal theory of transistors can conveniently 
be used to demonstrate the method. For a p—n—p transistor the governing 
equation for the flow of minority carriers in the base region can be written (1) 

2p p Pn l op (1) 
L* D, ct 
where p = the density of the minority carriers (holes) in the base region, 
p, = the equilibrium density of the holes in the base region, 
- the diffusion length ,(D,7,,), where D, is the relevant diffusion 
constant and 7,, is the volume lifetime. 


4 


p- 


(Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 2 (1956)] 
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For the purpose of calculating the emitter and collector currents we 
shall assume that in Fig. | the resistivity of the n-type base region is very 
much higher than the resistivity of the p-type regions and that therefore 
the potential difference across each junction is constant over the area of 
that junction. It follows that the emitter and collector currents flow 
normally to the surface of the junction. 

\ssuming abrupt boundaries at the emitter and collector junctions we 
have, for the hole density in the base region at the emitter, p, = p,, <'*? 


and at the collector p, p, <e*T, where } 


’ and J, denote respectively the 
potentials across the emitter and collector, e is the hole charge, 7’ is the 
temperature in degrees Kelvin, and k is Boltzmann’s constant. 

The recombination loss at the open surface of the base region leads to 


the boundary condition 


where n is the outward normal and s is the surface recombination velocity. 
In this paper the material of the base region will be considered to be 
homogeneous; but it is possible by the methods to be described to study 
the effects of arbitrarily positioned imperfections and variations in 7, and 
s over the surface and volume. 
The total emitter current is equal to 


I cD, | Kas I.nds,, (3) 


Ch 


where n is the positive outward normal (so that the partial derivative is 
egative for emitter current entering the base region) and S, is the area of 
the emitter junction. 


Similarly, the total collector current is equal to 


* an . 
] eD I dS, I.nds,, (4) 
2 Clt J 
where S. denotes the area of the collector junction. 


The total current lost due to surface recombination over the open surface 


the base re¢ o10nN is 
/ es P—Pp ) dSo | I.n dq, (5) 
ere S, denotes the area of the open surface. 
‘he continuity equation fo! the base region is 
) ), l 
puPe_-¥.1. (6) 


> f 


Hence the integral of the divergence of the hole current I taken over 
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the entire volume v of the base region is 


I, (p—p,,) dv—e | a dv | V.I dv. (7 
Tp . . C . 
The divergence theorem states that 
| V.Idv = | I.ndS,+ | I.ndS,+ | I.nd8,, (8) 
so that I. I+L- Ris (9) 


A relaxation solution of (1) (with the boundary conditions specified 
above) will show values of the hole density p at each mesh point. The 
integrands of (3), (4), (5) and (7) are all in terms of p and its derivatives 
multiplied by their respective constants. Relaxation nets can be caused 
to display at each mesh point on the surface the current lost by recombina- 
tion in each element of surface at each mesh point, and the divergence at 
each mesh point in the volume, merely by multiplying the relevant p values 
by the appropriate constants. 

veference to the governing equation (1) shows that the quantities shown 
in (9) will be functions of frequency. 

As usual in transistor analysis we use a quantity a which is a function of 
frequency. « can be defined by the relationship 


] -eD, | (ep/en) dS, 
x(f) E _ aa (10) 
_ D, | (eép/éen) dS, 
The current gain factor (2) is then 
 @ : (11 


‘ l—a 
for grounded emitter connexions of the transistor. 


3. Relaxation solutions for small signal conditions as a function of 

frequency 

To solve the governing equation (1) by relaxation we shall assume small- 
signal conditions. We substitute p’+-dp’e' for p in (1), where 6 < | and 
p’ isa complex quantity which we denote by a+ ib. The governing equation 
can then be split into a steady state expression and a small signal a.c. 
expression. The former expression is discarded and the small signal a.c. 
expression in turn splits into real and imaginary parts which, on cancelling 
out de and putting p, = 0, can be written 


: ' 
V2a = r3 (a—wz7, b) for the real part | 
? (12) 
” l ‘ : ; 
V°b I (6+w7,,a) for the imaginary part | 
#4 


where w denotes the angular frequency. 
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The real and imaginary parts of the emitter boundary condition become 


by a similar process respectively 


a, Pp lekT b, 0). (13) 
, The real and imaginary parts of the collector boundary condition become 
respectively (since V7 is large and negative during operation of the transistor) 
As 0, b, 0. (14) 
The real and imaginary parts of the open surface boundary condition 
? become respectively 
ca S cb s a 
a, b. (15) 
on D, on D, 
To calculate the frequency characteristics which are described later in 
this paper the equilibrium density p,, of the holes was taken to be zero (as 
» above). This was done because at zero frequency the a.c. solution is then 
identical with the steady state solution. This simplification involves no 
significant loss of accuracy of the steady state solution as far as our present 
purposes are concerned. 
The small signal a.c. forms of the integral expressions (3), (4), (5) and 
7) are found similarly to be 
r ca 
I, eD, | dS, 
ch 
| « 
4 cb 
I eD dS, 
| J Ch 
| 
F ca . 
| c. eD ds. 
J cn 
eo ee 
/ eD ds, . 
J cn . (16) 
} 
F. % adv CW b dv 7; t QQ, 
I] F bdv CW | a dv T, + Q, 
I, = es | adS, 
i es | bdS, 





where the suffixes 7 and 7 respectively indicate real and imaginary parts. 
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Using (9), (10), and (16) we can show that for the small signal case we 
have respectively for the real and imaginary parts of « 


Cf) wx der — lar + He) + Tel bet Ier+ Ie) ) 
™ Leet Le (17) 
vf) — Ler boi Heri + Hoi) Toi Ter Tore + Tor) | 
i B+, 
Then the phase angle of « is 
De Lu¢ — Loris +103) — Lei Lop — Leer +1. 
HD) = tan LL i 


The collector current J, does not appear explicitly in these expressions, 
This is an advantage in computing « because numerical differentiation is 
then required only for evaluating the emitter current /,,.+-i/ 


ei* 


The magni- 
tude of the small signal current gain factor is 
\ 2 


- | fotttad—a¥\? (a 
I / | X,( a a es ee = : lf 
. hoes? 2) (a ay)? + 3 | — 


Preliminary calculations showed that if a transistor is reasonably well 
designed (with the emitter area smaller than the collector area and with a 
reasonably small spacing between them) and if the material is as good as 
that normally used in manufacture, the hole density p is negligible except 
between the junctions and in their neighbourhood. There is then no hole 
density around the periphery of a square germanium wafer (Fig. 1). The 
emitter and collector areas used were symmetrical around an axis through 
the centre of the wafer and therefore radial coordinates (7, z) could be used. 
It was only necessary to compute the solution in one 7,z plane. This is 
equivalent to omitting the corners of the square germanium wafer, but no 
appreciable error is involved in doing this because the hole density in these 
corners is negligible. 

To solve a differential equation by relaxation this must first be converted 
into finite difference form. This process (3), and the process of solution by 
relaxation (4, 5) have been clearly explained in the literature and need not 
be repeated here. 


The finite difference form of the real and imaginary parts of (12) are 


h h h*)\ h? 
a,,+a_,,4 1 | 5) ; ( >,)" - (4 | za) +A, r ory bo | 
oat ” . “p 
h ji h? h? 
b,. 1 6 7 ( = i Ju { (1 >} - (4 | 73} >0 | A, p23 O70 
2) 29 4‘, f 4 


(20) 
where the A’s represent Fox’s difference terms (3, 5) and the notation is 
that used by him. The interval of tabulation is denoted by h. 
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\t the open surface boundary, we have (for the real and imaginary parts) 


h h h?— 2hs' h? 
2a_; (1 5) ar T (1 5, )@ Ir (+4 LD Ja +A; = — 1307») 
_ - “p p “p 


h h h? = 2hs h? 
Dy) l b,.+ {1 __- (4 |... )o LA = wT, a 
! | >) ; | >,) , a’ a," lC(<iCérK Sh 
(21) 


where n represents the outward positive normal to the open surfaces. These 
surfaces were arranged for convenience to lie along the r-axis. The collector 
and emitter boundary conditions are given by (13) and (14). 

Equations (20) and (21) can be regarded as referring to two separate 
relaxation nets for the two functions a and b, one of which is real and the 
other imaginary, which are inter-connected by the right-hand sides of these 
equations. The two nets are relaxed simultaneously. The boundary con- 
litions are such that when wr,, is equated to zero the imaginary solution 
vanishes leaving only the real solution. The process converges rapidly to 
1 solution unless wr, is very large indeed compared with L7,, when con- 


vergence becomes unsatisfactory. 


4. Examples of computed results 

For the purpose of studying the problems of design of transistors a 
number of differing shapes of the base region with various materials were 
computed so as to study the effects of these different shapes and material 
constants. The solutions of (12) were usually computed using three signifi- 
cant figures and one guarding figure. This represents a maximum error 
considerably less than is ordinarily attained in experimental measurements. 

Fig. 2 shows a relaxation representation of a typical fused impurity 
transistor. The boundary shapes of the emitter and collector can be 
represented for any shape by the usual relaxation artifices. The shapes 
shown were chosen in accordance with the experimental evidence at the 
time these calculations were performed. 

The figures in parenthesis represent current lost due to surface recom- 
ination in each radial element of surface at the position indicated. They 
must be multiplied by 10-® to represent amperes/em.? The figures not in 
parenthesis are the divergence current with respect to each radial element 
f volume; they must be multiplied by 10-4 to represent amperes/cm.? 
lhe original solutions used intervals of tabulation which were much smaller 
than those indicated in this figure and used graded meshes where necessary. 
lhe nets therefore covered a large area and it is not practicable to reproduce 
them in full in the space available. Only a few of the mesh points are shown 


n Fig 2 
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Table 1 shows the collected results of a calculation on the transistor of 
Fig. 2. Fig. 3 shows the corresponding frequency characteristic. 

It has been shown by R. L. Pritchard (2) that « and current gain for 
grounded emitter connexions as calculated analytically from the well- 
known one-dimensional solution of (1) are almost independent of the 
volume recombination time 7, at extreme frequencies. Table 1 shows 
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Fic. 3. A plot of the magnitude of the grounded emitter current gain against 


frequency for the transistor calculations shown in Table 1. 


that if the frequency in the three-dimensional case is made sufficiently 
high the terms Q, and Q, which contain w in (16) outweigh the terms 
7, and T', and also the surface recombination current /,,, and J,,, in the 
calculation of « and current gain. At extreme frequencies « and the 


i 


current gain then depend little upon either of the constants 7,, and s. 
This effect does not occur until the frequency is so high that the current 
gain nears unity. This extreme frequency much exceeds the cut-off fre- 
quency if that is defined as the frequency at which the current gain has 
fallen by more than 3 dB. One-dimensional and three-dimensional calcu- 
lations of the gain tend at such extreme frequencies to differ only by the 
lowering of gain in the three-dimensional case due to the spread of the 
hole flow away from the space between the junctions into the base. This 
spread is, of course, inherent in an actual fused impurity transistor. 

[t is found by studying the real and imaginary nets for each frequency 
that as the frequency increases to extreme values over the range of Table | 
there is a very slight decrease in the spread of the function p in the base. 
This decrease does not occur to an important extent within the range of 
frequencies of the transistor exemplified above. 


The effect of changing the geometry of a transistor was studied by 
: = ‘ . 


carrying out the integrations of equation (16) over different parts of the 
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base region for different geometries. The end to be sought in shaping the 


junctions is that the flow of holes from the emitter be confined as far as 


possible to the emitter surface immediately opposite the collector. During 
the process of calculating various cases it became obvious that not only 
does it follow that the collector should have a greater area than the emitter 
and be close to it, but that the area of the edge of the emitter which is 
exposed to the base region should be minimized as far as possible. If the 
material constants of, for example, the transistor shown in Fig. 2 and 
Table 1 are improved from s 5,000 em./see. to s 90 em./sec. and from 
7 10 microseconds to 7, 100 microseconds it was found, by a com- 
parison of the two solutions for zero frequency, that the spread of the 
p function away from the junctions into the base region was greater for 
the improved material than for the less effective material. It follows that 
the shape of the base region is a more important factor in the attainment of 
maximum possible gain when high quality base material is used than when 
lower quality material is used. 

The effect of increasing the overall size of a transistor does not leave the 
original distribution of the p function unchanged. In a typical instance an 
increase of linear dimensions of a solution by 5-906 times, leaving s and 7,, 
unaltered, caused the p function at zero frequency to extend all the way 
from the junctions to the edge of the wafer, whereas the spread for the 
original size was no more than comparable to that shown in Fig. 2; « fell 
from 0-9269 to 0-4806 and the magnitude of the current gain from 12-7 
to 0-925. The same changes in « aad current gain would have been found if 
the transistor had remained the same size, but 7,, had been reduced from 
100 microseconds to 2-87 microseconds and s increased from 90 cm./sec. 
to 531-5 em./seec. The importance of geometrical shaping and material 
constants becomes greater as size and power handling capacity increases. 

It has been found possible to invert the calculations. The problem then 
is not to find « and current gain when the geometry and material con- 
stants s and 7, are specified, but, for example, to find the value of s needed 
with a given geometry and a given value of 7,, to produce prescribed values 
of « and current gain. In the typical case of the geometry of Fig. 2 a 
current gain of 85-8, which was within an allowable error of 7 per cent. 
of the prescribed figure, was found in this way to correspond to s = 401 
em./sec. when 7, was 100 microseconds. The figure for s was arrived at in 
two approximations. 

It is possible to caleulate the saturated current which flows to the 
collector when it is at a potential J, greater than the saturation value and 
ly, is zero. This current is due to thermal excitation of the material. 


In this case p, cannot be equated to zero, but must be given a value 
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corresponding to the material and the temperature. Table 2 shows a 
typical calculation for the transistor of Table 1. The resistivity p of the 
materialwas taken to be 2ohms/em. which corresponds to room temperature, 


TABLE 2 


Zero emitter potential. Zero frequency. 
Collector potential more than saturated value 





Emitter current L, 1-166 pA 
Total net volume recombination and 

spontaneous generation current | £& 0-281 pA 
Total net surface recombination and 

generation current no. 1-197 pA 
Collector saturation current a 2-644 pA 





5. A comparison of methods 

Because analytical solutions of the transistor equations are only possible 
for a limited number of quite unreal boundary shapes, numerical analysis 
appears to be the only method of quantitatively analysing transistors with 
respect to the boundary shapes used in engineering practice. 

Because analytical solutions are so limited in scope an attempt has 
nevertheless been made by Moore and Pankove (6) to solve equation (1) 
for realistic shapes by measuring the flow of electric current through a sheet 
of conductive paper cut to the shape of a cross-section of the base region. 
Resistors were attached to the edges to simulate the loss of current by 
surface recombination. However, the current flow in a conductive sheet 
satisfies not the relevant governing equation (1), but the Laplacian Vp = 0. 
Because the term (1/D,)ép/ét is neglected, the paper model car simulate 
only the zero frequency case, and, because the term (p—>p,,)/7,, is ignored, 
the loss due to volume recombination is neglected resulting in an error in the 
p function and in the surface recombination integral /,,. The state of the 
technique of transistor manufacture did not permit Moore and Pankove to 
make any controlled comparison between measurements on an actual 
transistor and the results from a paper model. The present author has 
used the methods of the present paper to calculate at zero frequency 
some common shapes of transistor for typical ranges of s from 90 cm./sec. 
to 5,000 em./sec., and of 7, from 10 microseconds to 100 microseconds. 


Integration showed that the numerical value of the volume recombination 
current /, was never less than about a quarter of the total surface re- 
combination current /,. and was sometimes greater than /,,. Because « 
approaches unity and is equal at zero frequency to 1—(J 
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umerically negative) it appears that the paper model is greatly in error 


3 an inevitable consequence of neglecting J, altogether. 
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ON THE DYNAMICS OF PHASE GROWTH? 


By P. L. CHAMBRE 
(University of California, Berkeley, California, U.S.A.) 


| Received 23 June 1955} 


SUMMARY 
The growth of a solid, starting from negligible initial dimensions with either a 
plane, cylindrical, or spherical boundary, into a surrounding supercooled fluid is 
discussed. The densities of the solid and fluid are assumed to be unequal. As a 
consequence a convective motion occurs in the fluid having the characteristics of a 
source or sink flow because of the fact that a unit mass of fluid occupies, on solidi- 
fication, a volume differing from the volume originally occupied. The effects of this 


motion on the rate of growth of the solid phase have been analysed. 


1. Introduction 

THE study of the growth of a new phase in an unlimited uniform medium 
has been the subject of a number of investigations. One of the better known 
examples is the case of a plane ice-sheet advancing into supercooled water. 
Its rate of growth is essentially controlled by the conduction of the latent 
heat of solidification away from the boundary separating the two phases. 
A physically restricted form of this problem was originally solved by 
Neumann, and Ingersoll and Zobel (1) have dealt with certain extensions. 
Mathematically comparable solutions for cylindrical and spherical boun- 
daries were first published by Rieck and have been recently the subject of 
an extensive discussion by Frank (2). 

All the investigations mentioned above suffer from the physical limita- 
tion that the densities of both phases comprising the system are assumed 
equal. Hence, if one or both of the phases are fluids convection effects can 
be neglected. 

It is the purpose of this paper to analyse the dynamics of the growth in 
a two-phase system of unequal densities. The effects of convection must 
then be included and become important if the densities differ appreciably. 
The discussion to be given will deal quite generally with the growth of a 
solid, with either a plane, cylindrical, or spherical boundary, starting from 
initially negligible dimensions into a surrounding supercooled fluid. The 
convection in the latter modifies the advance of the solid boundary since 
it affects the ‘conduction of the latent heat away from the front accom- 
panying the phase change. Incidentally the results may equally well be 
applied to the diffusion of material if the phase growth depends on the 
latter. 

+ This work was performed under the sponsorship of the Office of Naval Research. 


(Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 2 (1956)] 
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Convective motion in the surrounding medium will arise if a unit mass 
f fluid occupies on solidification a volume differing from that originally 
ecupied. It will be assumed that this convective motion obeys the incom- 
pressible Navier-Stokes equations in a field of constant pressure. The latter 
sumption appears reasonable if one considers that the pressure difference 
ross either a cylindrical or spherical surface of separation of the two 
phases varies inversely with the radius of curvature of the surface. Hence 
the pressure difference ceases to be of importance once the solid has grown 
to a sufficient size. During the very early stages of the phenomenon the 
formation of a single nucleus of growth is believed to be due to fluctuations 
in the supercooled phase, Volmer (3). The macroscopic analysis does not 
:pply to this initial régime but becomes valid once the nucleus has grown 
sufficiently. The pressure difference is altogether absent in the case of a 
plane surfuce. The non-linear inertia terms are retained in the equations of 
motion since they become of importance for very large density differences. 
They must be considered in any case during the earlier stages of growth 
since large convective velocities arise in the fluid close to the solidification 


iront. 


2. Derivation of the equations 

The growing phase is assumed to be a solid of infinite thermal conduc- 
tivity so that it remains throughout its history at the constant temperature 
T. of its surface. Hence no transport equations are needed for the solid. 
Taking account of the above mentioned assumptions, the equations of 


motion and energy of the fluid take the following forms: 





U; UU V| Une t K{- U, “) (1) 
r rj 
rout.—"|r.+* 7], (2) 
a r 








where u(r,¢) and 7'(r,t) are the velocity and temperature at a point r at time 
!; and o, the Prandtl number, is defined by o = v/a where rv is the kinematic 
viscosity and a the thermal! diffusivity of the fluid. These parameters are 
assumed to be constants. A in equations (1) and (2) takes on the values 
), 1, 2 respectively for the plane, cylindrical, or spherical solidification 
front and suffixes denote derivatives. 

The boundary conditions for the convective velocity and temperature 


of the fluid at infinity are 


(00, t) Q, (3) 


T'(oo,t) = T 


x 


(constant). (4) 
Q 
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In addition the following condition must hold at the (as yet unknown) 
solidification front r = R(t), 
T(R,t) = T. (constant). (5) 
The principle of the conservation of mass yields a condition on this phase 
boundary for the convective velocity. The solidification front moves with 
a velocity dR/dt. Hence the mass solidified per unit surface area per unit 
time is p,dR/dt, where p, is the density of the solid. In the fluid there 
exists a convective motion which at the surface of separation of the two 
phases has a velocity u(R,t). Thus the net velocity is [d R/dt—u(R, t)| and 
this causes a mass flow per unit surface area per unit time of amount 
p,|d R/dt—u(R, t)| where p, is the density of the fluid. But 


dR dR 
“ - Rt = 2 p) 
A dt “\ ] a 
or u(R,t) = —e ak (6) 
where em 24. 


It is to be noted that the convection induced by the difference in density 
of the two phases ceases when the densities are equal. Furthermore, the 
convective flow will be directed towards or away from the solidification 
front according as «€ is positive (p, > p,) or negative. The value of « lies 
between —1 < «€ < ©, but a negative value implies that the density of 
the growing solid is less than that of the fluid phase and this is rare for 
physically realizable transformations. 

Finally, the thermal boundary condition requires that the heat liberated 
in the solid per unit area per unit time, be equal to the amount of heat 
transported through unit area of the solid phase boundary per unit time, 


so that 


= —A, TR, t)+ pi Q, T, — 
( 


where A,, C,, L are respectively the thermal conductivity, specific heat, and 


p(L-+-C, T.) 


) a u(R, 0}. 


heat of solidification at an appropriate reference temperature of the fluid, 
C, is the specific heat of the solid. Using (6) this simplifies to 


RD one | gtk (7) 
a dt 
where  - - 47. E |, a a ’ 
O C, Ci pi 


The equations (1) to (7) are the equations of the problem to be solved, 


subject to R(t) being specified. 
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3. A particular solution of the problem 

The above equations reduce to the conventional case of no convection 
ife, and hence w, are set equal to zero. A particular solution in this instance 
is obtained by prescribing a position—time relationship for the solid phase 
boundary of the form R(t) oc vt. It will be seen in what follows that a 
solution is possible in the presence of convection if one writes similarly 


R(t) = 2Bv(vt), (8) 


where 8 is a positive number to be determined. 
Similarity and invariant considerations of equations (1) to (8) suggest 


that the velocity and temperature have the following forms: 
/v 
u ‘4 — 4 (n), (9) 
T = Tn), (10) 


were 7) 


Del)’ (11) 
Hence both the convective velocity in the fluid and the velocity of the 
solidification front, which ave large during the early stages of the growth, 
vary inversely as the square root of the time. 

Substituting (9) into (1) reduces the partial differential equation in r 
and ¢ to the following ordinary differential equation for f(7) 


bf’+K(— f at) — ff’ +(f+nf’) = 0 (12) 


where the primes denote differentiation with respect to . A single inte- 
gration and use of the boundary condition (3), together with the assump- 
tion that f and f’ vanish at infinity, yields a Bernoulli equation 


2n)f < (13) 


\ 7 


which can be integrated in finite terms to give 
| (¢ s* sk) ds . (14) 


The remaining constant of integration A is evaluated from the boundary 
condition (6) which prescribes the convective velocity at the surface of 
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separation of the two phases. Since the position of this boundary is given 
by equation (8) one finds that 


uU <b /*. for 7 = B,t > 0, (15} 


which yields for the convective velocity distribution 


7 


Ifv\ e-™ | on, {2 es Sy Pree ie 
u = i : 7] — eBK+1er : + ¢BK+1eR* (e-* /e*) ” . (16) 
vis 7) 77 


For the different solidification front geometries one obtains, after some 
calculations, the following results which can be expressed in terms of the 
tabulated error function, erf (x), and the exponential integral, — Ei(—2), 


Plane front (K = 0) 


; 'o B? 
u —2 (2 ) . n| Tides 2h oo al a q (17) 
A \(2/ v7) + «Be? erf(m)—erf (8) || 








I/v\ e-?? «B2eP 
tae” (=) 7 CO CoEECLA ater es eat ” 
Spherical front (K 2) 
Uu 
, [(v\e* 93,62) > 3@(2 fe e-*) oe |" 
J (*) ; aha «B%e (= 2B ~| - 2[erf(n) erf(8)]}) } 
(19) 


These expressions hold for 7 > B, t > 0. 

We next calculate the temperature distribution; 7’ is given by equations 
(10) and (11) and together with the general velocity distribution (16) reduce 
the energy equation (2) to the following ordinary differential equation 
e-7" 


1" +07") 2n+ 


oT Na ok 


se) [ay 





9 ) ty s* 1 4 
BR =a eBK+1¢B? = : ds| +3 | — 0, 
Nir 


(20) 
The particular solution of this equation which satisfies the condition (4) at 
infinity is given by 
0 t 
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The boundary condition for the temperature gradient (7) determines the 
constant of integration A with the result that 


T_T 9 \20 2 9 2 ff e-8 oe dt 
2o(e+ 1 | - BK+1eop" e-P] — 4 eBK+1eh? — | —. ds —. 
0 Vir J | Nor vx J sk * 

7 
, (22) 


It is from this equation that the parameter f is calculated by imposing 
the remaining boundary condition (5) which takes the form 7(f) = T,. 
Equation (22) reveals that 8 depends upon the prescribed parameters 
T,—T..)/0, o, «, and K, and a complete discussion is obviously a problem 
of considerable complexity. Once 8 has been found the position of the 
solidification front, the velocity, and temperature distributions in the 
fluid are all completely determined by equations (8), (16), and (22), and 
the problem is solved. One can readily show that for the physically 

0, the velocity and temperature distributions 


- B,t> 0. 


tally if one neglects the non-linear inertia terms in the equation of motion, 


realizable range — 1 € 
given above are well defined bounded functions for 7 Inciden- 
which is justified for small density differences and large time, equations 
11) again lead to solutions. They are of simpler form than those 


4) to 


given abov e 


4. Discussion of the results for a fluid with Prandtl number 


unity 


The following discussion deals with the special case when the thermal 
diffusivity and kinematic viscosity of the fluid are equal. This is approxi- 
mately true for a gas. For a value of a 1 equation (22) can be integrated 


explicitly to ¥ ield 
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The remaining integrations can be carried out in terms of tabulated 


functions with the following results. 
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Since the case of the plane solidification front illustrates the essential 
properties of these solutions the subsequent discussion is restricted to 
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Fic. 2. The dimensionless velocity distribution for B o = 1. 
equations (17) and (24). The value of the parameter f in these equations 
is determined from 
T. 2(e+ 1)Be?| 1 —erf(B) i 
(a prescribed constant) = P if saat B)| (27) 
6 (2/ vz) + «Be*| 1 —erf(B)| 
One observes that solutions for positive 8 and —1 < « < © exist only for 
1 restricted range of (7.—7T.)/@. Since this parameter increases mono- 
) tonically with f for a fixed e one obtains the following asymptotic expression 
for a very large growth rate 
as 7 2(e+1) 


= » B>ti. (28) 
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This shows that for « = 0 one can, at most, prescribe 7|—T_, = 0. This 


amounts to the physical limitation that the condition of supercooling of 
the surrounding medium is just equal to that required for adiabatic freezing. 
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Fic. 3. The dimensionless temperature distribution for b=@ l. 


If « becomes positive and very large (7'—T,,)/@ cannot be assigned a value 

greater than two. On the other hand, for the opposite extreme, as ¢ tends 

towards minus one, a physically sensible solution becomes impossible. 
Fig. 1 shows a few curves of 8 plotted against (7.—7_.)/@ for different 


values of the parameter «. In the case « = } the density of the surrounding 


fluid is two-thirds the density of the solid phase, while for « = —} the 
density of the solid phase is two-thirds the density of the surrounding 
fluid. Similarly the density ratios for the values « = 20 and « = —2 


21 


supplement each other. In Fig. 2 the absolute magnitude of the dimension- 
less velocity distribution u(7,t)/u(B,t) is plotted against « for the specific 
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value of the growth rate parameter 8 = 1. For a large positive «, which 
corresponds to a strong flow directed towards the solidification front, 
the velocity profile is quite steep in contrast to an outward flow which 
corresponds to a negative e«. It should be kept in mind that the magnitude 
of the velocity u(1,t)/2,(v/7t) which forms the reference scale for these 
curves varies considerably with ¢ as shown in the insert of Fig. 2. Finally 
the dimensionless temperature distributions (7'— T,,)/(7,— T,,) are exhibited 
in Fig. 3 for a value of 8 = 1 and for different values of the parameter e. 
The effect of the convection is clearly discernible. The inward flow steepens 
the temperature profile in contrast to the outward flow. The maximum 





reference temperature |(7.— 7))/@|g_, and its dependence on « is shown in 


R 
yo 


the insert. 
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A VARIATIONAL PRINCIPLE FOR THE 
CONDUCTION OF HEAT 
By LL. G. CHAMBERS (University College of North Wales) 
[Received 21 April 1955] 


SUMMARY 


A variational principle is indicated for the equation of conduction of heat. 


1. Introduction 

A LARGE number of variational principles exist in mechanics and electro- 
magnetism, but there appear to be no such principles for the equation of 
heat conduction. However, recent work by Herivel (1) dealing with a 
dissipation function in which the variation is in terms of the generalized 
velocities only (the generalized coordinates being left undisturbed) indi- 
cates an appropriate line of attack, and it will be shown that a variational 
principle exists in heat conduction in which only the time derivative of 
the temperature is varied, the temperature itself remaining unchanged. 


2. Analysis 
Let V be a volume and let S be the bounding surface, 
6 = the temperature, 
t = the time, 
X = 06/ct = the rate of increase of temperature, 
H,, = the normal heat flow per unit area across S, 
q = rate of production of heat per unit volume within V’. 
K = heat conductivity of medium within V, 
p = density of medium within V, 
8s = specific heat of medium within V. 
[t is assumed throughout that there is no motion of the medium. 
Consider the following quantity and its variation with respect to X. 


r 2\2.. §«68(, Fo encl a ff Oo. fe 
Pu | as yar at | K(voy" dV — | qd | H, — as 
; J } 4 





x ct 

. dj (2.1) 
1 | psX? dV — K(V0).(WX) dV | qX dV H, X d8. 

V V r’ S (2.2) 


The functions @ and X are, from the point of view of variations, independent 


in the sense that in variational mechanics the quantities q,, 4, are indepen- 
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dent. The variation which follows is performed with respect to X only, 


and 
aF = [{ psX8X dV — | K(W6).(WSX) dV — [{ q8X dV + [ H,dX dS. 
r 8 (2.3) 
Using Kelvin’s generalization of Green’s theorem (2), 
SF [psX —q+V.(KV6)| 5X dV 4 | |, gs | 5x dS, (2.4) 
on 
j S 
whence it follows that if 
4 
H, Ks over the boundary, (2.5) 
Ch 
o6 = ee ‘ 
and ps ps.X q—V.(KV0@) within the medium, (2.6) 


5F vanishes. That these equations are satisfied follows from the fact that 
2.5) is the equation of heat conduction across a surface (3) and (2.6) is 
the equation of heat conduction within a medium where heat is being 
generated (4). 

It remains to be proved that the converse also holds, that is if 5 F vanishes, 
then the heat conduction equations follow as a consequence. ‘irst con- 
sider a variation 6X which, though arbitrary, vanishes on the boundary. 
Then ; -( 2 , ‘oe _ 

SF ps — o+'¥ (ey oe dV, (2.7) 


| C 


: 
and since 6X is arbitrary the quantity in brackets must vanish, giving the 
equation of heat conduction within a medium. If 6X is now arbitrary 
sr = | la,—K “\sxas, (2.8) 
z= on| 
s 
the volume integral being zero as its integrand has been shown to vanish 
identically. From this equation it follows that, as 6X is arbitrary, the 
vanishing of 5F implies the vanishing of the quantity {H,—K(e6/én)} 
which is equivalent to the equation of conduction of heat across a surface. 
The second varia‘ ‘on of F can be calculated easily and is clearly given by 
6°F =} | ps(dX)?dV > 0 (2.9) 
; 
whence it follows that F is a minimum for small variations in X (= 06/ét) 
about the true value. 
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AN APPROXIMATE METHOD FOR DETERMINING THE 
TEMPERATURES REACHED IN STEADY MOTION 
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SUMMARY 


rhe problem is resolved into two steps, one involving heat generation and trans- 


port, and the other conduction. A numerical and graphical method for solving the 
problem is given. As an illustrative example, the isothermals for cold extrusion 
through a smooth square die with a 66-7 per cent. reduction are found. The results 


are compared with the (analytical) solution obtained when conduction is neglected. 


1. Introduction 

THE plastic deformation of materials is characterized by irrecoverable 
strains. The work of deformation appears largely in the form of heat energy 
and indeed some of the early investigations into the nature of heat were 
prompted by an observation of the large amounts of heat generated in the 
boring of cannon (Rumford, 1). Farren and Taylor (2) measured the heat 
evolved as a percentage of the work expended in the tensile deformation 
of various metals. They found that for steel 86} per cent. of the work 
appeared as heat energy, 90}—92 per cent. for copper, 92—93 per cent. for 
polycrystalline aluminium, and 95-954 per cent. for aluminium single 
crystals. It is clear from these data that the temperature rise in a heavy 
deformation must be considerable. For many practical purposes this 
heating effect is unimportant, constituting at worst a nuisance. However, 
in other applications the modification of the material properties and the 
mechanics of the distortion may be serious. Thus the defect known as 
‘hot-shortness’ resulting in the circumferential cracking of the surface of 
extruded materials, appears to be attributable to the rise in temperature 
due to plastic working. In this paper the problem of evaluating the local 
temperatures existing in plastically deforming materials is considered. 


2. General considerations 

The general problem to be examined is that of time-dependent heat flow 
in an incompressible moving medium with heat generation in the medium. 
The velocity at any point will be determined from the solution of the 
problem of plastic flow and in general will vary from point to point. 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 2 (1956)] 
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(Conduction across boundaries separating different media must be con- 
sidered. Several analytical solutions exist to idealized problems of this 
class (see Carslaw and Jaeger, 3), but these are for simple boundary con- 
ditions. In this paper we shall be concerned more with practical problems, 
involving in many cases variables which have been determined numerically 
or graphically. The approach adopted here is therefore numerical. Only 
problems of steady-state plastic flow under plane-strain conditions with 
two-dimensional heat flow are specifically studied here. 

It will be assumed that a fraction «(0 < a < 1) ofthe energy of distortion 
is converted into heat energy. Modification of the plastic-flow problem due 
to the effects of temperature on the material properties is not taken into 
account. Let 2 and y be two orthogonal fixed directions. The differential 


equation to be satisfied is 


(2.1) 


- = U- v- : ae, 
Ox= = oy*| Ox cy Spc ot 


e706 a6) o6 cb aw 06 
where @ is the temperature, wu and v the velocity components in the x and y 
directions respectively, zw the rate of expenditure of plastic work per unit 
volume, and ¢ the time: «, J, p, and c are the thermal diffusivity, the 
mechanical equivalent of heat, the density, and the specific heat respec- 
tively. The rate of work w at a point is expressible in terms of the yield 
shear stress k and the velocity components u,, v, in the directions 2,, y, 


tangential to the curves of zero rate of extension at the point, viz. 


w = ky X(' = =) (: 


OY, Ox, 


to 
bo 
— 


u, v, and vw are therefore known functions of 2 and y. 

Problems governed by equation (2.1) are ‘one-point’ or ‘marching’ 
problems in ¢. Conditions given at t = 0 enable the problem to be solved 
step by step for given increments in t. Suppose the solution to a problem 
has been carried forward to some stage. Consider firstly a relaxation 
solution. Since equation (2.1) is of odd order in ¢ it cannot be directly 
converted to a difference equation for a relaxation attack. A change of 
variable must be made to make the equation of even order (i.e. second order) 
in the new variable, and an additional boundary condition nominated 
which does not conflict with the existing ones. Since there are three inde- 
perdent variables, the relaxation process must be carried out on a three- 
dimensional net. Finally, a numerical differentiation must be carried out 
to regain ¢ as one of the variables. The reader is referred to a paper by 
Allen (4) for a detailed discussion of this process. Since the first derivatives 
of @ with respect to x and y have variable coefficients, the relaxation pattern 
will vary from point to point. It therefore seems clear that a relaxation 
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attack would be extremely lengthy if not impracticable. The fact that 
equation (2.1) does not have constant coefficients is a serious practical 
disadvantage using any method, since it involves different numerical 
coefficients in the difference equation for each point of space-time. Ij 
became clear in the course of the present investigation that if a workable 
method were to be evolved, a considerable simplification of the problem 
would be necessary. It was therefore decided to separate the problem into 
two parts. 

Consider a small displacement increment of the plastically-deforming 
material. Simultaneously heat is being generated, material (and hence 
heat) transported, and conduction is occurring. The approximation to be 
used here consists of regarding the deformation as proceeding by a series 
of infinitely rapid jerks. Thus all the heat generation and transportation 
is regarded as occurring instantaneously followed by an interval in which 
conduction takes place as for a stationary medium. The length of this 
static period is the interval which would be required in the actual process 
for the material to traverse the distance travelled in the jerk. It would 
appear that by making the intervals sufficiently short the solution could 
be made as accurate as desired, except perhaps in the immediate neighbour- 
hood of a discontinuity line. 

There is no difficulty in principle in computing the heat generation and 
transportation. The problem is therefore to solve the equation 

, ev19 — (2.3) 
ot 

with @ a known function of x and y initially. Equation (2.3) can be put 
into a very convenient difference form by replacing both the time and 
space derivatives by their central difference relations. This form unfortu- 
nately leads to a rapid and uncontrollable build-up of rounding errors 
(Crank and Nicolson, 5). Hartree (6) gives a discussion of the problem and 
reviews methods which have been used to solve it. Some of these are quite 
straightforward but involve considerable labour. For the purposes of the 
present paper and in view of the approximations already made, it was felt 
that the work entailed in using the methods outlined by Hartree was not 
justifiable if not actually prohibitive. It was decided therefore to use the 

crude approximation to the time derivative 
et ie - [8, .(¢-+8t)—8, ,(t)], (2.4) 
where 6¢ is the interval in ¢ and 6, ,(t) is the temperature at the point (j,k) 
in x and y at time ¢. This expression does not lead to a build-up of rounding 
errors provided the intervals are suitably chosen, although the error of 
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approximation is much larger than that involved in the central difference 
expression. The difference equation incorporating (2.4) has been solved 
for various problems for which analytical solutions are known, and has 
been found to give results which are satisfactory for the purpose in hand. 


3. Details of the present method 
Problems of plane plastic strain are governed by equations which are 
hyperbolic in both stresses and velocities. The displacement or velocity 
field can therefore be represented in terms of ‘characteristics’ across which 
discontinuities in the tangential velocity component may occur. The 
velocity field is usually obtained either numerically (Hill, 7a) or graphically 
(Prager, 8). In certain cases it may be given analytically. The characteristic 
separating the rigid material from the plastically deforming material is 
commonly a velocity discontinuity. Thus, in general, the solution to the 
plastic problem will consist of regions in which no distortion is occurring, 
regions in which the velocity varies continuously and boundaries across 
which velocity discontinuities occur. It is remarked in passing that the 
velocity discontinuity separating rigid from plastically deforming material 
is only one particular type. At rigid constraints (e.g. a platen or die) slip 
will commonly occur and there will be an associated coefficient of friction. 
Consider a time interval 4¢. The motion of all points in the interval is deter- 
minate. In particular the initial location of all elements which are at 
arbitrary grid-points at the termination of the interval can be found. The 
amount of work done and hence the heat generated in the elements must 
also be found. If the element undergoes a continuous distortion, the rate 
of work is given by the formula 
CU, , Og Ch od 
1 f+ U, —— — Up ? 


O82 C84 C8, O8g 


wy = | (3.1) 


where s, and sg are elements of the characteristics with tangential velocity 
components uw, and vg respectively, and ¢ is the counterclockwise angle 
made by the «a-characteristics with a fixed direction. The reader is referred 
to Hill (9) for a detailed discussion of the plastic-working aspects of this 
subject. If the element passes through a region in which the rate of working 
varies rapidly, it may be necessary to evaluate the total work done in the 
interval d¢ in a series of steps. The plastic work expended on an element 
crossing a velocity discontinuity is performed instantaneously and can be 
found as the product of the shear stress k and the shear strain. The shear 
strain is equal to the magnitude of the velocity discontinuity divided by 
the normal velocity component at the discontinuity. At a constraint with 
friction, the work done can be found from the frictional force and the 
distance travelled in the interval. In general an element will receive heat 
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by a combination of these mechanisms. A ‘transport’ chart can therefore rule wou 

be constructed as follows. Divide the (x,y) or physical plane into a suitable | traverse 

network. As indicated above the initial location of, and the work expended | numeric 

on the elements arriving at the grid-points at the termination of an interval \ single st 

dt can be found. Assuming that no conduction occurs and that a fraction y 


steps if 


~~ 


of the work is converted into heat, a temperature rise can be associated } js not di 
with each grid-point. Thus each point in the net will have attached to it | ture (e. 
a vector, giving the initial position of the element arriving there in the flow vie 


interval 5¢, and a number giving the temperature rise occurring in its 


passage. > 4, Exa 

The heat conduction problem can be solved approximately as follows. ; Ag ar 
Using the approximation (2.4) for the time derivative and the central This is 
difference formula cent. re 

V6; x(t) Fal jsra(l +O; p(t) +9; pa r(t) +0; 4 -1(t) 40; x (0)], 
where h is the interval in x and y, (2.3) can be put in the form f 
P KT 

Oi e(t-+7) = 8; ,(t)+ 73 Disa a(t) +951 c(t) +8} pa r(t) +; 4-1(t)— 46; 2(t)], 

(3.2) 
where 7 is the time interval. 

Thus the value of @ at (j,/) at time (¢+-7) can be found from the solution 
at time ¢. 

Suppose the value of ¢ is known everywhere at time ¢ = 0. Consider a ; 
specific grid-point. The vector on the transport chart gives the initial ? 
location of the element arriving there in an interval 7, and by reference to 
the known solution, its initial temperature is known. Hence the tempera- 
ture of the element arriving at the grid-point at time t = 7 can be found by of the 
summing the (known) initial temperature and the temperature rise (whichis > Fig. : 


recorded on the chart). The conduction problem is now solved using equa- 
tion (3.2) for an interval of duration +. The repetition of the sequence 

(i) transport+ heat generation, 

(ii) static conduction 
enables an approximate solution of the problem to be found. Caution must 
be exercised in the choice of the time and space intervals. Normally velocity 
discontinuities and associated relatively high local temperatures will occur. 
Suppose the space intervals have been conveniently chosen. If the time 
interval is too short, the importance of the discontinuity will be magnified. 
All points arriving at the discontinuity will be credited with the heat 
generated (instantaneously) at it. This will imply that a zone of width 
equal to the space interval has attained the temperature associated with 
the discontinuity, a result which would be quite erroneous. A safe working 
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fore | rule would appear to be to choose the time interval such that the distance 


ible | traversed by the slowest element is of the order of the space interval. The 

ded | numerical solution of the conduction problem need not be achieved in a 

rval \ single step of duration r. It may be desirable to split it into two or more 

na | steps if the gradients are large. The incorporation of boundary conditions 
) 


is not discussed in detail here since it is adequately dealt with in the litera- 
O it ture (e.g. Allen, 4). For a discussion of boundary conditions from the heat 
the flow viewpoint, the reader is referred to Carslaw and Jaeger (3). 


4, Example 











™ | \s an illustration of the method, a problem in extrusion was considered. 
sae This is plane-strain extrusion through a smooth square die with a 66-7 per 
cent. reduction. The dimensions selected are shown in Fig. 1. The speed 
} hamber wa » 
2 y wah Va 
Scm DIE 
mn N 
nfec4 | Zero heat loss 
i BILLET R . a 
a , ; 
> Ycm 15 cm/sec. 
d 
: Le M s e a ene e e ee eae 
LO 
Fic. | 
1 
Vy {the ram is 5em./sec. The plastic solution to this problem is due to Hill (9). 
is Fig. 2 shows the slip-line field from which the velocities may be deduced. 
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In this particular example simple analytical expressions for the velocity 
components can be given. For the conduction problem it was decided to 
regard the billet as initially cold and with a high thermal diffusivity. The 
chamber and die were taken to have a low thermal diffusivity. These 
conditions would be met in the cold extrusion of copper through a steel 
die. It is appreciated that these conditions do not correspond to normal 
commercial practice for copper. They were chosen both to explore the 
practical application of the method, for which purpose initial conditions 
are immaterial, and to examine a pioblem in which conduction is a signifi- 
cant effect. Perfect contact was assumed to occur between the billet and 
the chamber and die, and no heat to be lost from the extruded material. 
The physical constants of the billet chamber and die are given in Table 1. 


TABLE | 





p ¢ K 
Billet 10 | oO! bare) 


Chamber and die | 10 | or | ov1 


« was taken to be unity and kg/Jpc to be 22-5° C. corresponding to a value 
of k of approximately 10 kg./mm.? Intervals of 0-5 em. in x and y witha 
transport interval of 0-1 sec. were taken. This time interval was resolved 
into smaller intervals (down to 0-025 sec.) depending on the temperature 
gradients. 


5. Discussion of example 
The isothermals at 0-3, 0-6, 0-9, and 1-2 seconds after the commencement 
of the extrusion are shown in Figs. 3—6. Consider the slip-line field in Fig. 2. 
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city the remaining material remains rigid. The ‘streamlines’ are indicated on 
to | the figure. An element located at P in Fig. 2 at the commencement of 
The | extrusion takes nearly | second to reach O during which time it suffers no 
ese } strain. Hence its temperature changes only by conduction. This is verified 
tel | in Figs. 3, 4, and 5. The main heating effects occur instantaneously on 
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crossing ABC and OC. Again this fact is confirmed in the isothermal 
patterns. The most interesting practical question to be answered is how 
conduction modifies the temperatures which would be achieved if the 
diffusivity were zero. 

The calculation of the temperatures which would be reached in the 
absence of conduction is not difficult. It involves merely the determination 
of the trajectory traced by each element and the plastic work involved in 
deforming it. Fig. 7 shows the isothermals determined by the current 

at . ° : 
method when conduction occurs, and analytically in the absence of con- 
1] 


duction. It is interesting to note that the exit surface temperatures are 
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substantially similar. The significant observation is that higher te 


mpera- 
tures are reached on the centre line and at the 


die-chamber wall corner when 
conduction occurs. This fact could be predicted on gener 
the nature of the velocity discontinuities. 
temperature distributions in the extruded m 
differ markedly. 


al grounds from 
As would be expected, the 
aterial remote from the die 
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If no conduction occurred, steady-state conditions would be achieved 
after one second. A review of the results suggests that after about 1:3 
seconds the isothermals to the left of the section PP’ in Fig. 


7 are nearly 
static. The temperature distribution in the die is 


again what would be 
expected from a knowledge of the deformation made of the billet. 
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Ta- that the temperatures fall off very rapidly. The temperature distribution 
hen in the die is not, of course, a steady-state distribution. A check on the 
‘om overall accuracy of the solution can be found by crmputing the total heat 
the } |caving the work material under steady-state conditions and comparing 
die it with the work done by the ram. If no loss of heat occurred to the die, 


these quantities should be equal. It appears that the current method over- 
estimates the heat evolved by about 10 per cent. The smoothness of the 
isothermals is surprisingly good, as is the consistency of values from one 
interval to the next in regions where steady-state conditions had been 
achieved. 

It will be remarked that the highest temperature in the extruded material 
occurs after the exit discontinuity has been crossed, a situation which is not 
physically possible. This is due to the motion being resolved into two 


distinct phases. 


: 6. Remarks 
The above example shows that the method suggested is workable. The 
relevant query is therefore to what extent the approximations made are 
justifiable and the results meaningful. Firstly, the method is considered 
from a mathematical viewpoint. Two main approximations have been 
made. The first is in considering the material to proceed by a series of jerks 
instead of smoothly and the second in replacing the differential equation by 
» adifference equation. The former ceases to be an approximation either when 
the time intervals are made vanishingly smallt+ or when the diffusivity is 
zero. Itisin any event concrollable. The chief error will be in the neighbour- 
hood of discontinuities. The approximation used for the time derivative 
in the difference equation is not one which would normally be employed 
except as a matter of expediency. It should, however, be remarked that 
material is only in regions of high temperature gradient for a comparatively 
short time. Moreover, if the time intervals are chosen sufficiently small, the 
errors introduced may be made as small as is desired. From the purely 
mathematical standpoint and in the absence of practical considerations, 
the methods suggested or reviewed by Hartree (6) are much to be preferred. 
A remark by Richardson (10) in a discussion of ‘Standards of Neglect’ 
would apr.ear to be relevant: “An error of form which is negligible in a 
haystack would be disastrous in a lens. Thus negligibility involves both 


mathematics and purpose.’ It is suggested that for the purposes of the 


ly | This is not true for the points crossing velocity discontinuities since they traverse the 
mtinuity in zero time. Hence no matter how short the time interval is, the temperature 

e 4 ; 2 : . 
associated with the discontinuity will be (relatively) large. Replacing the discontinuity 


ble the errors to be controlled. 
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present paper, the method described is adequate and that a reasonably 
accurate picture of the temperatures achieved both qualitatively and 
quantitatively has been found. 


The work described has been carried out as part of the research pro- 
gramme of the Mechanical Engineering Research Board of the Department 
of Scientific and Industrial Research. The paper is published by the 
permission of the Director of Mechanical Engineering Research. 
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\{ MODIFICATION OF LAGRANGE’S EQUATIONS FOR 
SMALL OSCILLATIONS WHEN SOME NATURAL 
FREQUENCIES ARE HIGH 
By HAROLD JEFFREYS (St. John’s College, Cambridge) 
Received 31 August 1955] 


Tue problem considered arises in such problems as the correction of the 
period of a pendulum for flexure of the rod, where the speeds of the free 
vibrations corresponding to distortions are very great compared with that 
of the rigid pendulum. Various methods have been used (1, 2). It has 
recurred recently in the treatment of the variations of the earth’s axis of 
rotation, where the internal elastic vibrations again have high free speeds 
and the corresponding coordinates can be regarded as giving only a small 
correction to the movements considered. In such conditions the contribu- 
tion of such internal motions to the potential energy may be appreciable, 
but terms containing the corresponding velocities will be negligible. We 
therefore consider problems where some of the velocities are absent from 
the Lagrangian. (The opposite case, where some of the coordinates are 
absent, can of course be treated by Routh’s well-known method.) 
We take n degrees of freedom and write the Lagrangian function as 


L L(q;,4;>t), (1) 
where, fori = k+ 1 ton, L does not contain g;. Then we have by Lagrange’s 
equations al 

. 0 («=k+1 ton). (2) 
cq 
' L 
Consider L’ L lq, (s = k+l1 ton). (3) 
~ OFs 
Then 
, c DL c L ° l oL ° a * L 
éL oq; 4 — Og; 00+ |e 0 - > 
eq, 64; 2\e4, 04, 
4 CLs ° l ( 4a ( Jae N OL) 
oq, . oq,-4 (" oq.4 gf de- q,8 — (4) 
cq, C4, ” \e4. 4s 4, 


where r runs from 1 to k and s from k+1 to n. 
Now (2) suffice to determine g, as functions of q,, ¢,, t. Suppose the 
variations of g, to be such that (2) are always satisfied however g, may vary. 


(Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 2 (1956)] 
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A MODIFICATION OF LAGRANGE’S EQUATIONS 


i 
rm él . OL , oL 
Then -= 0, 5— = 0, —_=0 (5) 
4s ds CFs CEN’ 
since g, does not occur in L. Hence the sum in brackets is zero. Hence if ‘ 
L’ is expressed in terms of q,, g,, t by using (2) to eliminate 4,, } 
éL’ @éL ol’ eb } 
— = > so > (6) By 
Cd; Cd, qd; 079; 
déL’ él’ 
and = (7) 


dt 04, eq, 
We notice that if Z contains quadratic and linear terms in the q,, the ’ 
former type are absent from L’ and the latter halved. We have thus a 
convenient method of calculating a modified Lagrangian from which the Un 
internal coordinates have been eliminated. Correct results would also be 
found if the factor 4} was omitted from (3), but it would be necessary to 
evaluate quadratic instead of linear terms in the q,, which is usually much 
more troublesome. 


The 
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SUMMARY 





the , aie ' . ‘ , 
[The central difference forms of the Euler transformation 
us : 
| a . 
he l u* d 5» OU m43 
th | l (1+ a)* 
» | 
. = 4 a)st! : . XY (1l—a)( ay" 52 : 
+ > 10°°U mis T ee Cine 
y Zaitt april t Z, 21+a)8+? t 
. s 1 8s 1 
uch 
ind 
| ea (—@)} ' — (l—a)(—a)*® ., 
, “ ] a enti a)*° LZ, 2%1+a)*%r 
s=1 s=1 
developed by means of the lozenge diagram from the forward difference form. 
[he latter formula is intended for use with asymptotic series whose terms alternate 
n sign, and to this end a table of the coefficients of the central and mean central 
differences in this formula is given. 
> 1. Introduction 
} Tue Euler transformation applicable to a series 5 (—a)"v, whose terms 
how n 
n 0 
ipproximate to those of a geometric progression of ratio —a, employing 
forward differences of the sequence v,, is well known. It may be derived 
formally as follows. We write 
Ev, a A E—1, 
t mm j 
’ then 
l ] ] 
Ugp—av, a'vs one , car age ¥ =~ Vo 
lt+ak l+a 1+ ja/(1+a)}A 
] a i ow 
¥9- 7 Avy T } A*v, mn 
lta| l+a \l-+-a, } 
In practice, use of the transformation is frequently delayed until the 
mth term, and we then have 
1 ] 
> (—a)"v,, = ¥ (—a)"v, +(—a) Om (1) 
n=0 = ltak 
J ( a)” a a 8 
S (—a)re,,4 5 Atv,,,. (2) 
hows It+a 4 \l-a 
be 0 s 0 


We shall develop formulae alternative to (2). 


[Quart, Journ. Mech. and Applied Math.. Vol. IX, Pt. 2 (1956)] 
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2. Formulae employing various differences 


where Vv =1—F-1 
we have 
] ] ] ] 1—V 


| ee a SS v : painless } 
l+ak ™  1+a 1+{a/(1+a){V/(I—V)} ™~ 14+a1—f1 (i-+a)}v ™ 
1 —+ ] n ] { ] n | 

" - es me Vn ) — > Vy . 
lta ( ) p3 (; ‘a “m = Tha ‘ioe . ( | ) Ym 


To obtain central difference forms and also to throw light upon their 
application we shall invoke the use of the lozenge diagram technique 





introduced by Fraser in connexion with interpolation formulae (1), and 
described by Whittaker and Robinson in (2). From the equations 


Sn—1y, —__Sn-ly — Sry 
041 — 9" Wy = "Vy 44, 





yfz™ti-e — ysz” 8 -gftize-s, 
—d l P 
where y= —, — ; 5 = E-+A. 
l+a ] +a 
we derive 
fe*-S*—*v, +-9'2***- 0, —_ ysz” shn Wnty’ +1 yn 8"v, ‘ae (3) 


If we now construct the lozenge: 
ysizn—s 6"-ly 
PS ysgnti-s 
aie 8"y 
ystizn—s p+} 
A 


Son —8 ly) / 
y®2" 8 \ 6” Un+1 


and travel from left to right along either the upper or lower paths, equation 
(3) shows that the sums of the quantities encountered on the paths are equal. 
Knowing, as we do from (2), the coefficients attached to a line of forward 
differences, we may build up a scheme of such lozenges and obtain the 
lozenge diagram: 


< Y \53,, ateye |) costs 
22y 22y2 | O'Un—} git lystt | oe Uma, 
Um \ 52, — j A Nd | 55 2y 
“Un m~ a+ > m 
wil oe w8+1))8 < ns 
\ zy y“\8 Unsa eae §2s+1, 
zy ¢ m+4 - zy” a di J -_ e8t+1y8+2) . 
zy” } o*t m+1 ie Y o*5 Um+41" 


snide | 
a bic ee one ere ae a glad 
zy3 o% m+grre oY grt nee” 


Expanding the operator in (1) in terms of the backward differences of c., 


one 


We 
exam] 
Gauss 


viz. 
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We are now in a position to obtain a variety of summation formulae, two 
examples of which are those using the schemes of differences found in the 
Gauss forward and Gauss backward interpolation formulae respectively, 


viz. 
. 25 29) lL sada) | aay2S52o, _| 429,253, 
Um AV» 17 a Um+2 — “Um ZYOV,, +h zy fc) Um ~ Y ) Um+} 
oy ~8+1),8+1528+1), a 48+1,,8+2528+2,, 
eUm > o Yy ) Um+3 p 3 = y 0 Um 
s=0 s=0 
and 
2; aT) ~ay0) w2ap529, _| »29,253., 1 
OO a AE ae Zn TH ZYOUpn — 4 TEYOUn P2°Y"O°Um_4 1 + 
9) 1 w8+1,)8+152s+1), 1 WY 48+29,8+1528+2,, 
Um > Y fc) Um A > rad Yy ) Um: 


s=0 s=0 


[t must be noted that the numerical application of these various series 
entails their curtailment at some term, either from considerations of ulti- 
mate convergence, or a heuristically inferred limiting relationship between 
successive terms and the remainders resulting from the neglect of all subse- 
quent terms. If we decide to traverse one route in the lozenge diagram and 
terminate the series so derived at one particular point, we may select an 
alternative route terminating at the same point and obtain precisely the 
same sum, since all series thus obtained involve the differences explicitly, 
and we are merely dramatizing a system of algebraic identities. There is 
thus no advantage in using any of the alternative summation formulae of 
this kind in preference to a formula employing forward differences. 


3. Formulae employing mean differences 
We now obtain expressions equivalent to the mean of two such derived 
formulae. For example, the Bessel form of the Euler transformation is 


given by 





Um avy, 1] avUm 2 
l a a 52 
Ve, - OU, 24+ = POV, 43 4 
Ita ™ (1+a)? ™™* (1+a)8 — 
2 3 
(1—a)a? ,. a . 
-- 0°v,..1— - - £O"V,, a — 
X1+ayt ™* Gays om 
| a. = (—a)tt 
ie ; aa 5 289; J 
ita m ( qe mts S a i qj H? Um+4 
g=] 





(1—a)(—a)**! 52841, 
+ = 2( ltaee° , a ae (4) 


a 
s=1 
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and the Stirling form by 
, ‘ iis 
Um — Wy TUVUp 2 
l a ‘ a (l—a)., 
——. n5v,,— ’ 522 
lta ™ (] a ya m a a )3 ) m 
a ~ a*(1—a) . 
- pv, T 3/1, bv nn 
(1+-a) 2(1-+-a)? 
] ~ (—ayr ., — (l—a)(—a)8 ,, 
ees Unc > pe Hote, \ > a ~ te e., (d) 
s=] ‘ s=1 
where p = 4( E+ E-4), 
For a = 1, (4) and (5) become 
Um —Um 1 -Um 2 

ve 15 a Nes eS _ = IT? 

5Um— tT Umip ts HOU m+; — 35 HOU. + 795 HOU msi —--- (6) 
—_—" = ; ley — pte, 1-5 ht — 1 de - 
and Un- Um+1 Um+2 — 2 Um ¢ HOV), 16 40 Um 34 HOU, sia dl (7) 


It is of interest to note that (4) and (5) may be derived with the help of an 
elegant artifice due to J. C. P. Miller. We write 


] . " 
v, = > C,8*%v,, + ¥ D8*v,,.,, 
| a E m jn s ¢ Um — Ss m+l1 
| - A? \ 
le. - - (C,.+D+D,A y. (8 
(1 it a) aX 2, 8 s s ) 1 | A ) 


Now multiply both sides of (8) by (1+<A)*-!, expand the remaining 
denominators in ascending powers of A, and equate coefficients of A®*—! and 
A*s; D,_, and C,+-D, are then respectively the coefficients of A**-! and A* 
in the expansion of (1+-A)s-1/(l1+a+aA). Thus 


ne alli — (’ Nts) teed 3 =) } | =" 


La (1+-a)*s 
a)s+ 
and C+D, a 
‘ (1 a)?s+1 
To derive (4) we write 
Um (Lu 58 )Un ays Um 1 (hu T 55)U h> 
then the coefficient of .6*%v,,,, in (4) is 
s+1 
C.+ D, ( pos 
, (1 | a)*8+1 
and the coefficient of 5**+1,,,, in (4) is 
. l (l—a)(—a)s+! 
4(D,—C,) = D,—4(D,+C,) = ( —— (e = I, 2...) 


2 (1-+a)?s+2 


When 4 





thus th 
—f (1 


then t 


and t 


Agai 
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When s 0, we have by inspection of (8), 


C42 , D, = - te Bs 
l+a (1++-a)* 


thus the coefficient of v,, in (4) is 1/(1+a) and the coefficient of dv,,,, is 


al(lta)?. To derive (5) we write 


then the coefficient of 6*5v, in (5) is 


C+ D,+4D,.1 


Again by inspecting (8) we see that the coefficient of v,, in (5) is 1/(1+-a). 


TABLE 1 





ry 5* 5° 5° 3° 
6289 
5 . } 
O24) 3245 ( 420) > 6940 
} 3896 7420 
7I4! 14300 
Id25 





Formula (7) is of particular interest when applied to asymptotic series 
whose terms alternate in sign, for it is characteristic of these series that 
the moduli of successive terms decrease to a minimum and then increase 
without limit. In the neighbourhood of the minimum term, therefore, the 
odd order differences of the moduli experience a change of sign, and the 
mean differences of odd order on a horizontal line through, or near to, this 
minimum term, will be relatively small. (Table 1, which displays the 
moduli and differences thereof relating to such an asymptotic series, 
illustrates this property In consequence of this property the notional 


convergence of (7) will be rapid. 
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An example of the power of (7) is afforded by the computation of For g 


2 — ly f (ap 1 se ] ) ~ r jn a/(1- 
o(%) = 4/($x7) {Jo(x) cos(x— 4a) + Yo(x) sin(x—}x)} ~ > (—1)"u,, 


r=0 of the C 
(4r—3)?(4r—1)? 
where u, = “Br 2r(8a Up—s M, = I, 








for the value x = 6:25. The sum of the first four terms ‘» tern 
| in te 


1—0-0018-+-0-00007 35—0-00000 96050 





is 0-99826 38950. Table 1 shows 10!w, and its differences for r = 4 to 11. 
Applying the forward form of the Euler transformation from r = 4, 

we have 

P,(6:25) = 0-99826 38950-+ 10-19 130444-3498-+ 1306-4 


+430+ 197+ 48+ 54—29-+4...) | 
— ()-99826 57527. | rs 


P4 


on stopping at the term before the smallest term, the last two figures 


being uncertain. f 
Applying the same transformation from the smallest term, viz. from 
r = 6, we find 
P,(6-25) = 0-99826 52941 + 10-19f4276 34 370— 165-4 223 — 285}. 
Transforming again, — 165+ 223— 285 + —{82—14+ 0} 
giving finally P,(6-25) = 0-99826 57516. c 
Applying the Stirling form of the Euler transformation from r = 7, 
we have 
P,(6-25) = 0-99826 61492— 10-1{4281 —373+ 64—7-L...} 
0-99826 57520, 
on stopping at the term before the smallest term, which is correct to ten 
decimals. 


4. The systematic application of mean central difference formulae 

If we have succeeded in applying (7) to an asymptotic series of descending 
powers of z, manifestly we may apply the transformation for some 2’ > z, 
and use (7) systematically to tabulate the function given by the asymptotic 
series for all z’ > z. This process, however, is laborious, since it entails the 
computation of sequences v}, and their subsequent differencing. Now in 
(5) the quantity a*/(1+-a)** is maximum for a = 1, thus if the transforma- 
tion succeeds for a = 1, it will do so for any a < 1, since the run of terms 
will decrease to a limit more rapidly. Moreover, for varying values of a, 
(5) entails the computation of but one sequence v,, and the differences 
thereof. It is thus more economical to write 1/2’ = a(1/z) (a < 1), work 


with a new independent variable a, and use (5) systematically. 
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of For general values of a, (5) may be computed effectively as a polynomial 
ina/(1-+-a)*. Tables 2 and 3 are given to indicate the orders of magnitude 
a 5 
of the coefficients 
A.(a) = a*/(1+a)?s, Ba) = (1—a)a’/2(1+a)*8+3 
in terms of which (5) may be written: 
| Um QV» ee Unis 
oll, | 
l 
4. | , : 8 S2s—1 8 S289, 
| Um > ( I )8A (a) 1d Um T > (— 1) B,(a) ot m* 
l+a g=1 s=1 
TABLE 2 
) 
10A, 10A, 1074, 10°4, 10°4, | 104A, | 104A, 
cures are) fe 00 0°6250 0°15625 0*3906 0°977 0°244 | 0°61 os 
00 249307 ‘62154 15496 3863 “963 "240 ‘60 “15 
8 s0914 *6og6¢ 15053 "3717 =| *gI5 °227 56 14 
from 7 42215 "58668 "14210 "3442 | °834 ‘202 oe; = 
( 375 "54932 12875 *3017 °707 -166 *39 "09 
49353 10974 *2439 *542 ‘120 27 “06 
4164 08500 1735 "354 "O72 ef | "03 
05594 "0993 176 031 | ‘06 | ol 
"02679 °0372 "052 oo7 | “cI co 
C 00564 0047 004 000 00 fete) 
00000 0*00000 0°0000 0*000 0°000 o"00 0°00 
‘. ” re - 
TABLE 3 
10B, o* Bs 10°B, 10°B, 10°B; 10°B, | 10°B, 
\oleie): oO > 0o*00000 0*0000 0°000 0*000 0°00 
ten 5607 - 04078 “1017 "253 063 “16 
7174 08263 | +2065 “SIO "126 | “31 
719 "12538 | +3037 736 178 | -43 
16093 *3772 “884 -207 | °49 
lae 370 { "18290 "4004 "903 *201 45 
. 1821 *3717 “754 "155 “32 
ing oe < 59 5: » 
15060 *2673 °475 "O54 "I5 
708931 "1240 "172 "024 03 
tic 02309 ‘OIgI “O16 “Oo! “00 
0o*00000 0*0000 0"000 0*000 0°00 
the ——— —________—— —__— 
in - ms - mee 
Using Tables 1, 2, and 3, we may now compute P,(6-58807 8459) 
la . 
corresponding to a = 0-9). The untransformed part of the series gives 
ms 
a, l (0-9) 0-OO01L8—+ (0-9)2 0-00007 35 
eS 
k 10-195 (0-9)3 96050 — (0-9)4 26088-+- (0-9)> 12097 — (0-9)® 8551} 
‘4 


0-99843 39847. 
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The transformed part of the series (from m = 7) gives 


af a 
10 10(0.9)7! 8563 — 0-2493 « 1492 —0-00656 « 2960-4 


{1-9 
L_()-0621 « 1024-+-0-00164 « 3245—0-0155 «x 426—0-00041 « 6940 


0-00000 01997, 
The total is thus 0-99843 37850, which is correct to ten decimals. 


5. Conclusion 

The requirement of a large number of figures in the final result precludes 
the early curtailment of the series in equations (4)—(7), and since v, is in 
general undefined for n < 0 (so that central differences of v,, after a certain 
order do not exist), the use of these formulae is invalidated, and this defect 
does seriously limit the extent to which these transformations may be 
applied. (Here it may be possible to apply the forward form of the Euler 
transformation to obtain a second asymptotic series, and so on, as in (3) 
and (4).) When applicable, however, formulae (5) and (7) do lead to con- 
siderable economy of effort. 
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